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GEOMETRICAL OPTICS OF ANGULAR STRATIFIED MEDIA* 
BY 
R. SEDNEY 
Aberdeen Proving Ground, Maryland 


Abstract. Some aspects of the geometrical optics of angular stratified media are 
considered. A similarity property of the light paths enables one to draw some general 
conclusions about the behaviour of light rays traversing a cone. An application to the 
interferometric method of observing supersonic flows shows that a conical flow test is 
valid including the refraction effect. 

Introduction. By the methods of geometrical optics, the paths of rays of light can be 
found by integrating a system of ordinary differential equations. The integration 
problem reduces to a quadrature in the two special cases when the isotropic medium is 
(i) plane stratified and (ii) spherically stratified. In these cases the index of refraction, n, 
is the function of a single variable, in (i) a Cartesian coordinate and in (ii) the distance 
from a point. In these cases it is possible to derive a number of properties of the rays 
without even specifying the functional form of n, [1] and [2]. In particular, the light 
rays remain in a plane so that these two cases are essentially two-dimensional problems. 
In this paper we consider media called angular stratified media, with the property that 
n is a function of the polar angle with respect to some line. Here the two- and three- 
dimensional problems are essentially different; however they both have a “similarity 
property”’. This similarity property was called to the author’s attention by Dr. J. H. 
Giese. 

In neither the two- nor three-dimensional cases was it found that the integration 
could be reduced to a quadrature. But some general conclusions regarding the light 
paths have been obtained. In the two-dimensional case the problem can be reduced 
to a quadrature plus a first order differential equation. 

Consideration of angular stratified media arises naturally in the study of supersonic 
flow over a cone or a wedge (Taylor-Maccoll or Prandtl-Meyer flow). Since optical 
methods of observing flows are quite common, the study of angular stratified media is 
not entirely academic. The paper concludes with an application to the interferometric 
method of observing supersonic flow over a cone. 

Three-dimensional case. The path of a ray of light in a medium of index of refrac- 
tion n can be written 


d(n dr/ds)/ds = grad n, (1) 


where r: (x, y, z) is the position vector of a point on the ray and s is the are length. 
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Equations (1) are not independent; in addition 

(dr/ds)? = 1. (2) 
Suppose n is a homogeneous function of degree zero 

n(Ar) = n(x), 


i.e., 2 is constant on radial lines. Then grad n is a vector whose components are homo- 
geneous functions of degree — 1. Then it is easy to see that (1) and (2) are invariant 
under the similarity transformation 


r—NMr, 
(3) 
s — As. 
This is the similarity property of the rays. An angular stratified medium is one in which 
n is homogeneous and rotationally symmetric, i.e., 


n= n(x + y*)'”*/2], (4) 


where z is the axis of symmetry. 
For definiteness consider the following optical problem. A cone of half-angle ¢, 
vertex at the origin, and axis along z, 
2+y’ =2' tan’ ¢ 
is imbedded in a space of constant index of refraction N, (see Fig. 1). There may be a 
discontinuity of index across the conical surface. Inside the cone the index is of the 


form (4). 


ry 








Fia. 1. 
A parallel beam of light travelling in the positive x direction strikes the cone. The ray 
tracing and image problems are to be considered. In the latter the mapping of points 
in some object plane x < 0 onto some image plane x > 0 is to be studied. 
In the ray tracing problem (1) and (2) must be integrated within the conical region 
subject to initial conditions for the position and direction of the ray. The initial direction 
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is obtained by applying Snell’s law. Consider the points, in an object plane, z = x, < 0, 
along a line J, through the origin, Fig. 2. 

The rays through the points of |, strike the cone along a generator. Thus the initial 
positions of these rays are related by the similarity transformation r — Ar. The initial 
2S 


Y, P 








Fic. 2. Object plane. 


direction, from Snell’s law, is invariant under this transformation since this direction 
is determined by the direction before refraction and the normal to the cone surface. 
Therefore, because of the similarity property of the rays and the invariance of initial 
direction, the light rays are similar curves. If r, and r, are the position vectors of two 
rays which strike the cone along a generator, then r, = Ar, ; for the geometrical path 
length at corresponding points, s. = As, . Furthermore the optical path lengths [Z,] 


and [L,], where 
[L] = [ nds, 
1) 
of these two rays are related by 


[L2] = [Ly]. (5) 


Tracing the paths of all rays which strike along a generator it is seen that these 
rays lie on a surface. This surface can be generated by the motion of a radial line. In 
particular, all these rays terminate along a generator and therefore all leave the cone 
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Fic. 3. Image plane. 


in the same direction. We have the following picture. The plane sheet of rays defining 
l, in the object plane, after traversing the cone, is still a plane at some angle to the original 
plane. Thus in an image plane x = x; > 0 the line J, is mapped into a straight line 1; , in 


general, not through the origin, Fig. 3. 
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Under certain conditions, e.g.,n > N, this mapping will be singular over certain regions 
of the image plane. The Jacobian of the transformation will be zero along a certain 
curve, the envelope of the lines J; . Between this envelope and the line S no rays are 
imaged. This accounts for the ‘“‘shadow of the shock wave’’ in a shadowgraph of super- 
sonic flow over a cone. Also, forn > N, the tangential rays, i.e., when the line J, is S, after 
traversing the cone do not emerge; they are internally reflected. This will also be true 
of a small bundle of rays in the neighbourhood of the tangential rays. The last two 
statements concerning the singularity of the mapping and internal reflection are easily 
proved for the special case n = constant and have been verified by numerical integration 
of (1) and (2) for a few cases with variable n. 

There is another result that can be obtained without resort to numerical integration. 
A lower bound on the depth of penetration of a ray into the conical region can be obtained. 
In [2] it is shown that, for any rotationally symmetric medium, the equations of the light 
rays can be transformed into the equations for a two-dimensional problem provided 
that the index of refraction is replaced by m, 


so 


m = [n’ — (h’/p’)]' 
where p = x” + 7’ and h is a constant for a given ray. This constant can be written as 
h = xq — yp, 
where p and g are the optical direction cosines with respect to the z and y axis. For our 

parallel bundle of rays, evaluating h before the ray enters the cone 
h= —yN, 
where 7, is the y coordinate at the piercing point. To obtain a real value of m 
p >h'/n? > N*y3/n3, 
where 7, is the maximum value of n = n(p/z). For example, in supersonic flow over a 
cone 7 is a monotone decreasing function of p/z so that ny is the index evaluated at the 


surface of the cone. In the few numerical cases considered this lower bound was found 
to be quite close to the point of deepest penetration. 

Two-dimensional case. A two-dimensional problem is obtained as a special case of 
the three-dimensional problem by considering a parallel sheet of rays in the x — z plane. 
In this case introduce polar coordinate r, 6. Then for an angular stratified medium 
n = n(@). It is evident that the similarity property holds also in this case. Using 6 as 
the independent variable instead of arc length, the path of a light ray is described by 
the equation 

[nr’ /(r? + 7°)'?) — nr/(r? +7)'? = 0, (6) 


where ’ = d/dé@. Let w be the angle between a radial line and the tangent to the path. 


Then 


2\1/2 


cosy=r/(r*’ +r)”, 
sin yy = ar/(r? + 7°)”. 
Then (6) can be written 
y’ — (n’/n) cosy = +1, (7) 


r'/r = cot ¥. (8) 
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The integration problem reduces to solving (7) plus the quadrature (8). It does not 
appear possible to solve (7) without specifying n. (The integration problem for (7) 
would be simple for special functions n, e.g., n = ae’’.) 

An application. In the interferometric method of observing supersonic flow, one 
tries to relate the density in the flow field to the shift in a fringe pattern. Since this is a 
tedius data-handling problem, it is useful to devise methods of getting partial infor- 
mation from a photograph of the fringe shifts. One such method is a test for conical 
flow, [3], which is a result of showing that the fringe shift is a homogeneous function of 
degree one in y and z. In [3] this was shown under the assumption that the light rays 
are not refracted as they pass through the conical flow region. The same result can be 
obtained allowing refraction. 

In the usual experimental arrangement a lens is placed in the one beam of light that 
passes through the disturbed region. The lens is placed so that the image of the median 
plane, x = 0, is photographed. Straight lines through the origin of an object plane, l, , 
are mapped into straight lines through the origin of this image plane. 

The fringe shift is proportional to the difference in optical path lengths of two inter- 
fering rays, one of which passes through the disturbance, the other does not. It has 
already been shown that for a ray through an angular stratified medium the optical 
path length has the similarity property (5). It remains to show that the same holds for 
the optical path length of the interfering, undisturbed ray. The interfering ray can be 
obtained as follows. The parallel bundle of undeviated rays can be projected back as a 
virtual beam (the interferometer actually separates this beam from the deviated rays). 
Taking account of the lens, it is necessary to identify an undeviated virtual ray with a 
deviated one that appears to come from the same point in the median plane (see Fig. 4) 
which is a projection of the three-dimensional phenomena onto a plane. 


VIRTUAL RAY 





DEVIATED RAY 








ME DIAN PLANE LENS 
Fia. 4. 


The appropriate undeviated ray is the one that passes through Q. Lay off a spherical 
segment with center Q and radius QP, . By Fermat’s principle the optical path lengths 
from P, and P{ to the point at which interference takes place are the same and also the 
optical path lengths from the light source to Py) and P% are the same. (Actually the 











230 R. SEDNEY [Vol. XIV, No. 3 


interferometer is usually arranged so that there is a constant difference between these 
optical path lengths but this is of no consequence in our considerations.) The difference 
in optical path lengths between P,P, and P,P{ must be shown to have the similarity 
property. From (5), [PoP,] has this property. Since Pj and P{ are determined by the 
initial and terminal points of the ray P,P, , the coordinates of Pj and P{ also have the 
similarity property. Since the virtual beam passes through only a region of constant 
index N 
[PoPi] = N PoP; . 

If 6 = [P.P,] — [PP] then 6 has the similarity property. Finally the coordinates in the 
median plane (y,, , 2,,) also have the similarity property. Therefore, 

SOY m » AZm) = AS(Ym y Zm)- (9) 
This is the basis for the conical flow test of [3]. Note that it is essential in proving (9) 
that the lens be used to focus on the median plane. For any other plane (9) would no 
longer be true. 
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A SOLUTION OF THE HEAT TRANSFER EQUATION FOR LAMINAR 
FLOW BETWEEN PARALLEL PLATES* 


BY 
Ss. C. R. DENNIS AND G. POOTS 
The Queen’s University of Belfast 


Introduction. Consider a viscous incompressible fluid flowing from left to right 
between infinite parallel boundary walls. The walls have a temperature @, to the left 
of the origin and zero to the right of it. If the heat transmissivity of the walls is taken 
to be negligible compared with that of the fluid and heat dissipation within the fluid is 
neglected, then for Poiseuille flow the dimensionless reduced temperature 3 (= 0/0.) 
satisfies the equation 


a8 a2) _ oe 
(x2 + an’ saa aE (1) 
with 
. 1 cs 2 
u(n) = Full — 75), (2) 


where « is the thermometric conductivity, u the longitudinal velocity component, 
uo the mean velocity and 2D the distance between the plates. If we introduce the dimen- 
sionless variables x and y defined by ¢ = 3Pé Dz, » = D(2y/m — 1) then with the sub- 
stitution (2), Eq. (1) becomes 


229 40% 10 (,_ va, = 
* Ox oy lr Y~ +) ax’ «~~ 3x Pé 





and if the fluid to the left of the origin is assumed to have the same temperature as the 
plates, the boundary conditions are 


0=0 for z>0, y = Oor7n, (3b) 
3’=1 for «=0 and 8 —>0 as tr~>o~, for O<y<z. 


The parameter Pé is the Péclet number and may be quite large, even for moderately 
small Reynolds numbers, so that for most practical purposes the first term in (3a) 
may be omitted. In fact in any case in which this term is not omitted, it is probably not 
justifiable to neglect the heat dissipation term in formulating (1). This term involves 
no significant change in the present mode of solution, however, and will not be con- 
sidered further. It should also be noted that, with the above simplification in (3a), and 
taking x as a time coordinate, this equation assumes a particular form of the heat con- 
duction equation for one-dimensional heat fiow in a medium whose thermal properties 
vary with position (but not with time or temperature). 

This paper is not concerned with thermal properties, but with a mathematical method 
of solution of these problems. The solution of (3a), subject to the conditions (3b), has 
been obtained when a is small by Prins, Mulder and Schenk [1], who solved using power 


*Received June 21, 1955. 
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series the corresponding Sturm-Liouville equation obtained on separation of variables. 
The method appears to be laborious, however, and our object is to indicate a new approach 
in which formal Fourier analysis, of the type developed by Carslaw and Jaeger [2], 
may be extended to problems of this kind. A detailed solution is given in the case of 
(3a) and the method is then extended to the general case in which u(y) may be any 
function satisfying Dirichlet’s conditions. It is capable of extension to three-dimensional 
problems, and, in particular, to the problem of forced heat convection through a duct 
with rectangular cross-section. 

Solution of the heat transfer equation by Fourier analysis. From consideration of 
the boundary conditions (3b) the local temperature distribution is symmetrical about 
y = 7/2 and is zero at y = 0 or for x > 0, so that we may expand #(z, y) as the Fourier 


series, 


A ib 


2D r9(2) sin ny (n = 1,3, 5, ---). (4) 
=1 


Then in the usual manner, multiplying each side of (3a) by sin ny and integrating with 
respect to y from y = 0 to y = x we obtain 


a” o% —nv, = = | (y - v’) O8 sin ny dy (5) 
dx *® Jo X aw / Ox 
with 
Z 
v,(0) = af v,(o) = 0. (5a) 
Putting 
yy” 1 
y-L=15,4 Sd, cosy (¢ =2,4,6,---), 
T 2 q=2 
where 
b, = = b, = - 
then we find 
"(, — ¥) 98. zis) d@, 8s. __p dv, 6 
I (y T Neale ada (z * 2rn°/ dx T : (n° — p’)* dx 9), (0) 


where the summation extends only to odd values of p. 
Finally we arrive at the infinite set of simultaneous constant coefficient equations 


in v,(z) given by 
2 ‘ 9 wood , 
= - (8 $,)% Os Li - We, (n ~ p) — nV, = 0 (7) 
dz? 3r wn 


dx T i (n? — p’)* dx 
and their solution is required subject to the conditions (5a). This system will be satisfied by 


v,(z) = A, exp (Az) (8) 
if 
128\n 5~ __pA, _ 
f(a, nA, + 7 4G —p) (n ~ p) = (9) 











~ 


= 
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where 
ai 8 8 
f(A, n) =a’ - (5 os $.), —n (10) 
3r rn 
and } is obtained from the condition for non-vanishing A, . This condition may be 
written compactly as 
A(A) = | a;,(A) | = 0 (¢,j = 1,3, 5, ---), (11)* 
where 


, - 
a oa # i), au = f0, i). 


a;; = 4;, = 3 
17 7 (2 aa 

The diagonal elements of (11) are very much larger, by virtue of the term 7’, than the 

non-diagonal elements [i.e. a;;(A) >> a@;,;(A)] thus as a first approximation to the solution 

of (11) we assume that 


) 


La;(\) | = T] anm(A) (12) 

m=1 
giving approximately the mth root X,, to be the solution of f(A,, , m) = 0. This equation 
has two real roots of opposite sign but only the negative one is acceptable on account 


of (5a). Corresponding to each X,, for m = 1, 3, 5 --~ there is an associated set of A<” 
forn = 1, 3, 5 --- and using the above approximation as an initial estimate of ,, a 


rapid process of successive approximation may be set up for calculating the correct 4,, 
and associated A,,”’ as follows. The mth equation of (9) is used to determine \,, . If we 
put A,, arbitrarily equal to unity in (9) we obtain for X,, the equation 


: 128m — = 
fn, m) + we * > int a AS” (m # p) = 0. (13) 
The remainder of (9), excluding n = m, may be written 
‘ wa) 128 = p (m) 128 mn 
1 ( , nA, = Aw = «i aes By atone 14 
fm sn)” + are. — pie Ae” (n ~ p) zt Ae Gat — ape (CA) 


and once X,, is given these may be used for determining the associated A,”’. The process 
is now as follows. Taking the initial estimate of \,, as the negative root of f(A,,, m) = 0 
we substitute this in (14). These equations are well conditioned because of the leading 
diagonal elements and very accurate estimates of the A{” may be obtained by neglecting 
the terms under the summation sign giving 


my _ _ 128 mn Fj . 
i = = Ke in? — np Lm 2). (15) 


From the estimates (15) the summations in (14) are computed and in this way we arrive 
at the A;”’ corresponding to the assumed \,, . We now return to (13) and re-calculate 
\,, and the process is repeated until sufficient accuracy is attained. In the numerical 





*In order to attach any meaning to (11) it must first be shown that the infinite determinant A(A) can 
be made convergent. Absolute convergence is readily achieved in this case by dividing each row by its 
corresponding leading diagonal element, provided that » is not such as to make any of these elements 
zero (see for example, Whittaker and Watson [3)). 
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examples considered it was found that if we aim at five decimal accuracy two complete 
iterations were sufficient. For example, in the case of large Pé, when the first term of 


(3a) may be neglected, the initial estimate of \43 = —32.2217. Using this value in (15) 
a set of A,” was obtained which were then improved using (14). The second estimate 
of A; = —32.1478 and on using (14) again only minor corrections occurred. The final 
accepted value of A; was found to be A; = —32.1472. 
It remains now to calculate the constants of integration as these have so far been 
chosen arbitrarily by setting A,” = 1. Thus, if we write 
2< i 
W(z,y) = =- @,, exp (A,,2)X,.(y), (16) 
T m=1 
where 


ro) 


X..(y) = D0 A,” sin py, 


p=1 
then using the boundary condition (3b) for v,(0) we have 
— 2 


Ax” Qn = 7 (n =I, 2, 5, e*s) (17) 
1 


- 


m= 


so that (17) provide an infinite set of simultaneous equations for the calculation of the 
@,, . For computational purposes however, these are unsatisfactory because of the slow 
convergence of the right hand side and it is better to proceed as follows. We make use 
of the variational principle that the quantities @,, must be so chosen that, when x = 0, 
the integral 


° © 2 
J(Q; ,@s,°°*) = [ (x= @,+4, — Il 2) dy 


/0 n=1 


shall be a minimum. Thus for m = 1, 3, 5, --+ we have 


oJ . - . 
—— = 2 | (x @,X, — I 2)X. dy = 0 
OQ m 0 n=1 

and on simplification of this integral we obtain the alternative set of simultaneous 
equations for the @,, given by 


ro) 


gS 24% o TF a,( As as”) (n = 1,3, 5, °°). (18) 


p=1 P p=1 
Although (18) are more complicated, the derivation of the @,, from them is in fact more 
accurate. In practice to quote the local temperature distribution in the regions of interest 
to within 0.1% it is necessary to obtain only the first four or five X,, since these occur 
in 3(z, y) as exp (A,,x). In the special case of the foregoing for which the term a’d°3/dx’ 
may be neglected in (3a), explicit formulae may be derived for the @,, on account of the 
orthogonality of the functions X,,(y). Thus 


® 2 
[ (y ~ xX.) =0 for m#n, 


ra 2 al 2 
a= Ff (v-“)x. av/ [ (v—“)xt@ ay (19) 


so that 
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and putting 
X,.(y) = >> AS” sin py, 
p= 


—— > > _— re Aim)? 1 ye 
Qn = a = 24 [. D ] _T 2 n? ’ 


p=1 P p= 


we obtain 





where 


co =) 


Lh” = LAMA = |p-a)—- DL VAPAP@=pt+a. (20) 
=1 q=l1 


p=1 q=l 


As a trial calculation, we have computed the first four eigenfunctions for the case 
treated by Prins, Mulder and Schenk [1] and the results are shown in Table 1. In this 


TABLE 1 




















| | 

m 1 3 | 5 | 7 
wn —2.82776 | —32.1472 | —93.4775 | —186.781 
Qn +1.8489 | 40.4860 |  +0.2407 +0. 1352 

p AS» | 4? | AS Fi 
+1.00000 | +0. 18615 | 10564 + .08356 
3 — 02118 | +1.00000 | + .42332 + .30054 
5 ~ 0014 | — 15617 |  +1.00000 + .66385 
7 | — 00020 | — .00429 | — .35303 +1.00000 
9 — .00006 | — 00171 | + .01157 — .62145 
11 — .00002 | — .00062 | — .00327 + .06758 
13 — .00001 | — .00027 | — .00132 — .00736 
15 — .00013 — .00042 — .00197 
17 | — 00007 | — 00035 — .00109 
19 — .00004 — .00020 — .00062 
2 | | — .00002 | — .00007 — .00038 
23 | — 00001 — .00004 — .00024 
25 — .00003 — .00016 
27 — .00002 — .00010 
29 — .00002 — .00007 
31 — .00005 
33 | — .00004 
35 | — .00003 











case the term in a’ is neglected, but it has been verified that, in such cases in which 
this term must be included, the calculations are more rapid and less A{” are required 
for each eigenfunction. In these cases, however, the dissipation term (£, 7) should 
be included in (1); the only effect of this is to introduce a known function of n on the 
right hand sides of (7). It should also be pointed out that, in the computational process, 
it was not found necessary to use the Rayleigh principle in estimating the eigenvalues, 
although this would undoubtedly give more accurate estimates at any stage. 
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Generalisation of the method. The method may readily be extended to the more 
general case of (1) where u(y) may be any function satisfying Dirichlet’s conditions. 
After appropriate transformation we should arrive at a modified (3a) with (16/z°) 
(y — y’/mr) replaced by some function g(y) which must now be expanded as a Fourier 


cosine series in (0, 7). 


i ‘ 
gy) =5 bo + >» b, cos qy (q = 1, 2, 3, ---). (21) 

- q=1 
The simultaneous differential equations can be formulated in terms of the coefficients, 
b, and the substitution (8) leads to the infinite determinant (11) where now 


Q;; = G;, = 4(b;.;, — 5);-¢)), a;; = f(r, 2) 
with 
f(r, 1) = a’d” — 4(bo — b2,) — @. 


The infinite determinant A(A) may be made convergent in the manner previously adopted 
[see footnote to Eq. (11)] provided that the sum of the off-diagonal terms a,; is absolutely 
convergent. If this is not the case, in fact provided only | a,; | is less than a positive 
constant, we may put A, = B,/n* in (8) and form the determinant for non-vanishing 
B,, , which can be made to converge as previously. 

It is probable that the determination of the eigenvalues and the solution of the simul- 
taneous equations may be performed in the manner described in the particular example, 
but the rapidity of convergence of the iterative procedures will of course depend upon 
the rate of decrease of the coefficients b, with increasing g. In any problem, approxima- 
tions to the eigenvalues \,, may be determined independently from the A,” by expanding 
A(A) about the top left hand corner and taking into account sufficient columns and rows. 
As some check on the numerical work, this has been done in the specific example, taking 
into account four columns and rows. These yielded approximations to the first four \,, 
given by \, = —2.82776, A; = —32.1482, 4; = —93.4968 and A; = —227.164. We do 
not consider other forms of the boundary conditions (3b) in this paper, except to point 
out that no significant change is involved by replacing 3(0, y) = 1 by 8(0, y) = a specified 
function of y. For example in the problem of flow between parallel plates, if the plates 
have respective temperatures +06, to the left of the origin, then we have #(0, y) = 
1 — (2y/x), and 3(z, y) may be represented by (4) with n = 2, 4,6 --- . It should also 
be pointed out that mathematically the method is basically equivalent to separation of 
variables and conversion of the resulting Sturm-Liouville equation to one of Hill’s type [4]. 
The mode of solution is not the same, however, and no restrictions, save those stated, 


need be placed on g(y). 
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ON THE INTEGRATION OF NON-LINEAR PARABOLIC EQUATIONS 
BY IMPLICIT DIFFERENCE METHODS* 
BY 
MILTON E. ROSE 
Office of Naval Research,! Washington, D. C. 


Abstract. The object of this paper is to investigate how solutions of mixed initial- 
boundary value problems for a certain class of non-linear parabolic equations may be 
obtained by solving suitable implicit difference equations on a rectangular lattice and 
taking the limit of such solutions as the mesh of the lattice tends to zero. 

We consider the non-linear’ parabolic differential equation 


Ou ( ou oe) 
al t,t,U (1) 


in the strip 0 < ¢ < T,0 < zx < L, with the initial condition 
—B'(x)u(x, 0) = fo(x) (1a) 
and the boundary conditions 


Ou 


a'(t) — (0, t) — B'(Hu(0, ) = f,(d, 


Ox (1 b) 
—a%(t) 2 u(L, t) — BOUL, ) = fuld, 


assuming that the solution u(z, ¢) is unique and exists with suitable regularity properties 
in the strip. 

Introducing a family of rectangular lattices 7, with mesh (h, k) and non-negative 
weights w, , w, , We associate with (1) a family of implicit difference equations 


w, Viu(x, t) + wV_u(z, t — k) (2) 


= F[zx, t, u(x, t), o, V(r, d) + w.Vau(2, t — k), Vila, dD) 


with corresponding difference equations approximating (la) and (1b) [see (1.1)]. 

In the first part of this paper we prove the pointwise convergence, as we let h, k — 0 
in such a way that k/h? = dis a fixed number, of the solutions u*(zx, ¢) of the difference 
equations (2) to u(x, t). We obtain estimates of the degree of convergence by methods 
suggested by arguments of Gévrey [7] and Laasonen [8] for linear parabolic differential 
equations and related to [4], [5], [6] for solving elliptic equations using a maximum 
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2For higher space dimensions our analysis would carry over with (1) replaced by 

Viu = F(t, 2, y, °°* 5 Us Us, Uz» Uy, °°*)- 
We will retain the form of our arguments suitable for extentions to higher dimensions. 
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principle. Our proof of convergence imposes a restriction on the value of \ [see (2.3)] 
which is stronger than the von Neumann criterion for linear equations. 

An iteration method is then discussed for solving the implicit difference equations 
(2); a convergence proof is given which provides both an existence theorem for the 
solution of (2) for a fixed value of h and estimates of the error at any stage of the iterations. 

1. Statement of the problem. Consider the non-linear parabolic operator 


ab , dd db 
Lio] = 2% — F(z, 4,4, % , %). 
le] = 3,3 1's Pr on at 
F(z, t, z, p, g) denotes a fixed continuous function of its variables for (z, ¢) in a region 
R in the (2, ¢)-plane and for all z, p, g. We assume that the partial derivative F’, , F, , F, 
exist, are continuous, and satisfy the inequalities 
O<a <F,<S A<@, 
IF, |<B<e, 
Ose <F,<C<o, 
where a’, c’, A, B, and C are fixed constants. 

Let R = R(T) be a domain in R bounded by the coordinate lines x = 0, ¢ = 0 and 
the lines x = L, ¢ = T; the closure of FR will be denoted by R’. The set composed of the 
segments B,(0 <x < L,t =0), B(x = 0,0 <t< T) and B,(z = L,0 <¢t < T) wil 
be denoted by B = B(T) and called the boundary of RP’. 

We define boundary operators A, , A, , Az by 


Aoly] = —B°(x)e(x, 0) on Bo , 
Ale] = a(t) 2 40, !) — B\()e(0, t) on B, 
2; 0 / 42 / y . 
A.[e] = —a’(d = g(L,t) — B()e(L, t) onB, ; 


here 6°, 8’ , 8° are continuous positive functions, a’, a are continuous non-negative 
functions on B, , B, , B, . Let fo , f: , fe be fixed functions defined on B, , B, , B, re- 
spectively. 

A mixed initial-boundary value problem ® may be formulated as follows: for fixed 
T, determine a function u(z, t) defined in R(T) with certain regularity properties satis- 
fying the equation 


Llu] = 0 inR 
and the initial and boundary conditions 
A,fu] = f; onB,, ¢=0,1,2. 


We will assume it to be known that this problem has at most one solution which exists 
with suitable regularity properties’ under appropriate regularity conditions on the 
operators L, A and on the initial and boundary data. 





‘Specifically, we assume Usrer , Uszt , Ut and lower order mixed partial derivatives exist and are 
continuous in R. 
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Before describing a method for approximating the solution of problem @ by differences 
it will be convenient to introduce certain notations. Let R,,, be a rectangular lattice 
covering FR’ given by lines 

x = mh, m=0,1,---,M, 

t = nk, n=0,1,---,N, 
where h = L/M,k = T/N. A net point P with coordinates mh, nk will by denoted by 
P,,., or, simply, (m, n). Let A, denote the set of points with coordinates (m, n) for 
0 <_m < M and let Af = A, + (0, n) + (M, n). A point of the set 


N 
Rix _— 2 A, 


n=1 

will be called an interior point of R,,, . A point of the set B,., = Ri. — R,.,isa boundary 
point of Ri, . Also, if w is a function defined on Ri, we will write w,,,, for w(Pa,.)- 
Furthermore, when F{,, is representative of a sequence of lattices obtained by letting 
h, k approach zero in such a way that the ratio 

k 

i? 
is fixed we write R,, = R,, By, = B, , ete. 

Next, let 


= 


Li lem x = w; Vem n + w2V em .n-1 _ F|z, ’ a ’ Omn » 01 V min + wW2V Pm.n-1 , VT @m.nly 
where 


1 
Viena _ h? (Pm+i.n + Pm-i.n — Minds 
ee eee (1.1) 
zPm.n = Dh Pm+ivn PYm-i.n)» ° 


1 
V nin - 5 Pmn es oe F 
Here w, > 0, w, > 0 and w, + w, = 1. Also, let 


Ax,ol@m.o} = —Biem.o ’ 0 < m < M, 
1 


an y 
Ax, 1[¢0.n] = h (Gi.n — ¥o.n) = Bivo.n ’ 0 < n < N, 


2 
An 
Ay olem.nl = h (Gn-1,.9 a €n.n) _ Bivn.n . 


The mixed initial-boundary value problem @, consists in finding a function u* defined 
on R, which satisfies the equation 
Lu} =0 mR, (1.2) 
and the initial and boundary conditions 


A,,;{u‘] = ft onB,.,, (1.3) 


A: . . 
where f; is a given function on B,,; . 
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In Sec. 4 we describe an effective procedure for computing u* and obtain an existence 
and uniqueness theorem for problem @, . In the next sections we derive sufficient con- 
ditions for solutions of @, to approximate the solution of ® uniformly in RP’. 

2. Estimates for linear difference operators. It will be sufficient to derive estimates 
for linear operators of the form 


xXr[em.n] = @1V Om n TT Wo V Omn-1 = fo a ane 


(2.1) 
— baal VPa.e + 2V Ha.n-1] — Cn.wPa.n » 
where 
<0 54.. $ 4.14.15 3.06 Shia SC. 
We assume h is so small, say h < hy , that 
Bh < Bh, < 26, 
where 0 < 6 < 1. Also, let \ = k/h’. 
Lemma 1. Let ¢ be defined on Rj with 
and xe] 20 mR, (2.2) 
A,.ile] > 0 onB,.; , t= 0,1, 2. 
Then, if 
wh .. (2.3) 


we have ov < Oin Ri. 


Proof. Let R; , B;,; denote the subsets of R, , B,,; respectively with points having 
t coordinates less than or equal to t, . Let M;[y] = max; max,,,, ¢, 7 = 0, 1, 2. We first 
show that 


maxr, ¢ < M\[¢]. (2.4) 


We proceed by induction and suppose that (2.4) is true for R, and By, , 1 <v < N, 
i = 0, 1, 2 and that there exists at least one point, say P,,,,,, , of A,., for which ¢,,,,.; = 
S > M;.'[¢]. From (2.1), (2.2) we obtain, since x,[¢] > 0, 


Om.n+10m,n+1 < ee | ee eer — Om.n+i) + ) ae oer On aaa) 
+ Lace Oues n a laan + le 1.nOm-i.n 
= | ee n+l 9 


where 


. h’ 
a 5 = WA ] + b,, nt+1l"o 9 


\ - 


Fa a. (Amr <— 2rAw>2) , 


/ 
lnc. = wd(I + Dane's). 








\- 
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By assumption, 
h | Oansi | < Bh < 26, 
so that 


: > wir = 6) Ps 0 


and Smnsicn + Lm-i.n = 2Aw, . From 2.3 we have alsoT,,,, > 0, Mnsi.n > Oand Pasiin + 
Preis + Pain = On.nvi > O. 


Suppose w, ~ 0. Then, by assumption, there exists at least one point P,,-,,.; of 
A,.; at which ¢,,:...,; = S and either (¢,-4:..41 — S) or (@m--1..+1 — S) or both are less 
than zero. If w, ¥ 0, 


, n Y 
max (Pm? +10 9 Pm’ .n ee < M;, x S. 
In any case, therefore, 
Gat 0a* acai — hen: n+1°5, 


which is impossible. For n = 1, i.e., Ri,, By,. (2.4) is obviously true and hence for 
Loan cc N. 


We next show that M;[¢] < 0. Suppose in fact ¢ had a positive maximum at a point 
Q of B,.; . If @ C By.o we would have 


Ax.ole(Q)] = —B°Q)¢e(Q) = o(Q) = 0, 


which is a contradiction if g(Q) > 0, since 6°(Q) > 0. If Q C B,., or B,.. and if Q’ is 
the h-neighbor of Q in R, , (2.4) shows that 


9(Q") — ¢(Q) < 0. 
In either case we would have 
0 < a(Q) — An. a,2[e(Q)] < —B"?(Q)e(Q), 


which is a contradiction if g(Q) > 0. Thus ¢(Q) < 0 for Q C B,.; , i = 0, 1, 2. This 
establishes the lemma. 

From this lemma and the fact that x, , A,,; are linear operators it follows that the 
only solution of the homogeneous problem is g = 0. The existence and uniqueness of a 
solution of the non-homogeneous problem follows directly from this fact. 


Lemma 2. If ¢', ¢’ are arbitrary functions on R{ with 
xal¢'] < —| xal¢"] 
and 
Ax. :[¢'] $=] Ay. s[¢"] I, 7=0,1,2, 
then 
le |<e¢'. 

The proof follows by letting ¢ = ¢’ — ¢’ and applying Lemma 1. 

3. Proof of convergence. If¢',y’ are arbitrary functions and if g = ¢' — ¢’, we have 


L,[¢'] = L,[¢"] _ xale], 
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‘where x, is given by (2.1) with 
Os,5 = Pas , ty , 26,4, ODs,7 Hei 5-1 , M4) 
and with b,,; = F,,c;; = F, evaluated at the same values of the arguments. Here 
24 = (rls tl — reader, 
Pi = (ri VAei J +0 — 1.) Veil}, 
G49 = (riiVileisd) +0 — 7.) Vile}, 
withO < 7;; <1. 


t.? 


h h h . ° 
Let v° = u* — u, where uw’, u are solutions of @, , ® respectively. Then 


xso"] = —L,[u). 


Now set 
(m, , Mz , M3 , M,, Ms) = MAX (| Urs |, | Use |, | Ure |, | Use |, | ere |) 
R’ 
and put 
mM, Ahm, 
A,(h) = 2 + 3 
and 
Bm 
A, = w.(m, + Ams) + —>-: 
From 
| L,{u] | = L,{u] 7 Llu] ly 


the usual Taylor expansion yields the estimate 
| Lalu] | < h?A,(h) + kA, . 


Similarly, if 


a* = max max | a’ |, M, = max max | u,, |, 
i Bi(T) i Bi(T) 
we find 
max | A, ,[u] — A,[u] | < ha*m,. (3.1) 


e ° A ° ° . + 
To obtain estimates for v” we introduce a comparison function 2* defined on R’ for 


which 


lA 
top 


Xr Q 
and 


A,[Q < —}. 





‘Explicit constructions of 2 are readily obtained following the construction in Bers [4]. 
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Clearly, 2 > 0 in R’. Let 
Q’ = max 2. 
Ry 
From (1.3), (3.1) we have 


max | A,,;[v"] | < (a*m,)h + max max | f; — fi |. 
i Bi(T) 


Let e > O be chosen so that 


€ > 2 max [(h?A,(h) + kA,), (a*moh + max | f — f”* |). 


B(T) 
Then 


xaleQ] < —5 < —| xa") | 


and 


< —| Ay; fo") |, ¢ = 0, 1, 2. 


Hence, from Lemma 2, we obtain 
h 
lv | < eM’, 
‘ e ° of . h 
which furnishes a uniform estimate for v’. 
Assuming for the moment the existence of solutions u* of ®, we may restate our 
results in the 


Let problem @ be approximated by problems @, in the sense that 


* - * 
max max | f; — fi | = Ja*O(h) if a* #0, 
i B;(T) } O(h’) if at a 0. 


Theorem. 


Then, under the conditions detailed above, in particular, if h, k — 0 in such a way that 


ional. (3.2) 


h? an 2w2 
the solutions u’ of ®, approximate the solution u of problem @ uniformly in R’(T); i.e., 
lu — wu | = a*O(h) + O(h’) + O(K). 


The condition on the mesh ratio \ given here differs from the condition for stability 
in the von Neumann sense for linear equations, where one would expect stability for 
arbitrary \ when w, < 4 (see [2], [3]). These facts raise the question whether our con- 
dition (3.2) may be weakened even in the non-linear case. 

R. Richtmyer has given® a convergence proof for the solution of the heat equation 
u, = Uz, by implicit difference equations of the type described here in which the von 
Neumann criterion is sufficient for his convergence arguments. From another viewpoint, 
condition (3.2) was imposed in order to obtain a maximum principle; F. John has shown 





5Seminar in Numerical Analysis, Institute of Math. Sciences, New York University, 1954. 
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11], at least for explicit difference methods, convergence without using assumptions of 
positively weighted coefficients. 

4. Construction of solutions of ©, by iterations. In this section we shall describe an 
effective iteration method for calculating solutions of ®, and obtain, at the same time, 
a constructive proof for the existence of solutions u* of @, . 

Let 


Gyii4(Z, or, S25 Vi) = olf: + S2 — 2Z) + w(Viars + Yi-n.; — 2Y,,) 


satel P| ’ ties 2,3 (oi - $2) + = (Via. = Y;-1 i), Fae = ‘|. 
Zn 2h y 


v 


Let ¢° satisfy the initial conditions 


ch 


A,ole ] = fo onB, 


and be otherwise given arbitrary on Rj . 
Let 0; be an ordering of the set A} . For p;,; C Af and, with 


Hi. |Z; t’, t?; Z*] = Z + 86,.(Z, t: , te 3 Y:.;) 


for 5 > 0, set 


¢; i = Hi. il). » (G)i+1.¢ » @i-1.5 5 Gi-1], (4-1) 
where 
(ints = Ginis 
if 7 + 1 precedes in the ordering 0; ; otherwise 
(G)int.s = Giar,s - 
At boundary points Q set 
-h 
ASQ) (4.2) 


n+1 a n 
D = COO ( ’) , — 
¢” (Q) Prag ¢)'(Q’) ate? 
Q’ being the interior h-neighbor of Q. 
Let ¢ = ¢° on A, in (4.1). We shall show that 


4 n 
lim ¢3,.1 > ¢i.1 - 


ei 
With ¢, , considered as initial values on Aj calculate ¢;,, from (4.1), (4.2). Continuing 
in this manner we will obtain a function ¢; ; defined in R, which, for P;.; C R, , satisfies 


y ° oa 
Gs, 5410@s.541 9 Pit+1,74+1 9 Pi-1,7+1 * ;.:) _— 0 


with 
A . 
A, slg] = fi on B,,; , a= 0, 1, 2, 
i.e. g;,; is a solution, by construction, of @,. 
Finally we shall prove that, given any e > 0, and for N = T'/k, there exist integers 
N,, M2, °*** , Ny depending on ¢ such that 


le. —~GiWiT)|<e @)CKh, 
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where the notation ¢7'; [g?/;"] indicates the function obtained on Aj} after n, iterations 
with “‘initial conditions” g7‘;2, on Aj_, . (We will write ¢7,; for ¢7,; [g;-,] where no con- 


fusion is apt to arise.) This permits the effective computation of the solution numerically. 
First of all, we have 


et — ets = (Stor - ort (FE) ia, - ory 


+ (FE) tas - OTL, 


the asterisk indicating that the partial derivatives are evaluated at intermediate points 


(4.3) 








a = —] 1 
hina ~ 1S, Ph 
0< 7, <1. 
At a boundary point Q (from 4.2) 
n+1 n ne ae mS 4 n , me n=l , 
°@Q) -—Q => + Bh [()"(Q’) — &)""(Q’)]. (4.4) 
Let 
m’,*' = max | ¢3";' — ¢i,; | 
ai 
and 
us’ = max |¢""'(Q) — ¢(Q) |. (4.5) 
(A’j-A4j) 
From (4.3) we have 
5 |x 8 |x | aH®. |* 
m:** < max |m? , uw] max {eet 4 | ee + | ie | \ (4.6) 
Now 
bx 2 
Oe et + i( 2 +ers + © F?), 
Oz k 
aH ** — ( h 's) 
at, = §w,| 1 9 >); 
ox 
a bu ( 1 +2 Ft). 
Of» 2 


( ‘hoose 


; h? -1 
5< (20, + h?A + MB) 


Then, since by assumption 1 — (h/2) F, > 0, all these derivatives are non-negative and 


2* he * 
(2 + sh + 2), ye =1- owes +“ rs) > > 0. 

















; 
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i 
Let ; 
' 
H 





p=1-— 6r°(A + k''B). 
Then (4.6) yields 
m’*" < pmax |m’*, yp"). 


From (4.4), however, we see that 


so that | 


or 
n+1 n 1 
m Spm;. 
Thus 
m:p" 
max | ¢777 — ¢3,; | < 7: (4.7) 
4a’; .~ p 
This implies the existence of a limit function ¢,,; , 
—— , ; . 
¢i.4 = lim ¢j,; , j= 1,2, -°> ,N, 


n--o 


which is, by construction, a solution of ®, . The proof that this solution is unique follows 
directly. 

Estimates similar to (4.7) hold for iterations based on a function w;_, taken as 
“initial values” on A/_, . Writing ¢;,; = ¢;,; [¢;-:], we have that 


| gis [ws] — gis | S | i [wi-i] — ¢.;[wj-1] | + | ess[wi-] — 


Also 
oH’. \* oH*.\* 
¢:.;[W;-1] — G5 = ( = ) (¢;,;[wy-1] — 9.5) + (“a ) (Gi+1,s[Wj-1] — Gi+1.5) 
\ dz Of) 
dH? ,\*, 
+ (24: (¢;-1,;[Wj-1] = Gi-1,5) 
Of» 


4. 


-{- 


9H’ ,\* ; ‘ 


OH; ig 
ay... (Wi41,4=1 —™ Pi+i,j-1) 


oH? + 
+ (4 (5-1, 7-1 —™ Yi-1,i )) 
é—1, 


the derivatives being evaluated at some intermediate point. 


Let 


U;[w;-.] = max | ¢.s[ws-s] =~ Pig | 
b’; 
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and 


V;[wy-1] = max | W;,;-1 — i.5-1 |. 
A’jn-s 


r) + 2 
(sy = fm. 20), 
of) . ( h :) 
(ee = bul — 5 FS), 


oH)" 7 ( h :) 
(fe = dw,| | + 5 FS 4 


U;[w;-1] < pU;[w;-1] + «V;[w;-,], 


Now 


Thus, from (4.8), 


where 
<A tas *y, 
oi: #8 a > hy, 
This may be written as 
U, [ws] $ = Vilws-l- (4.9) 


Now, let w;.; = ¢7'; where n; is some integer. From (4.7), 


ye ae 
Viejo S$; ox 5 Vilei-tl;, 
so that, with 
W? = max | ¢7,;[e¢520"] — ¢.;l¢75"'] (4.10) 
as 


we have 
W; < p'W;. (4.11) 
Thus, from (4.7), (4.9), (4.10) we obtain 


| eis lei')] — gs | S Wi + U lei’ 


i-l 


< p-Wi + 





(1 z pee Vile] j= 1, 2, hos ,N. 


We may assume V[go] = 0. Let 





K 0 
Tr; = max 3 V;[¢;- 
‘” 2 a-e 
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and assume 7, , 72, -** , My are so large that 


max W;|<T. 


2. 
j7=1 2,°**,N 


Let 
r = max (TI, , L,]. 
Then, 
¢7lejes')] — o.5 | S Te” + 0), sss FR 
and 
Gi.3— Gar] S Tp”. 


Hence, let « > 0 be given. Let n; be determined, for example, such that 


ess, 521 
with 
pe’ < ire 
Then 
T(p"' + p"'-") < say (2i aD ce 


and Tp” < e. 
This, then, furnishes a uniform approximation to the solution ¢;,; in R, , i.e., to the 


solution z' of @, 
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A METHOD FOR THE CONSTRUCTION OF GREEN’S FUNCTIONS* 


BY 
BRUNO A. BOLEY 


Columbia University 


Summary. A method is outlined for the determination of the Green’s function 
associated with any partial differential equation for arbitrary domains. The Green’s 
function is obtained as the solution of an integral equation. A method of solution of this 
equation is discussed which yields the Green’s function as the limit of an infinite sequence 
of functions. Convergence of this sequence is proved for the case of Helmholtz’ equation. 
An example from the theory of heat conduction is solved in detail. 


1. Introduction. The use of Green’s and related functions provides one of the most 
powerful tools for the solution of many boundary value problems of theoretical physics 
and applied mechanics [1, 2]. Its chief drawback lies in the difficulty of determining 
the Green’s function, for the particular partial differential equation concerned, in the 
case of domains with complicated geometrical configurations. For certain differential 
equations general methods are available; for instance, the Green’s function of the 
Helmholtz equation 


aa t jy? = Mes ule, W320 (1) 
can be derived for an arbitrary domain as the solution of an integral equation [1]. For 
boundary value problems not governed by Eq. (1), general methods are usually not 
available, and it would therefore be desirable to devise a procedure for the construction 
of Green’s and similar functions whose basic principles are independent of any specific 
partial differential equation. A procedure of this type is described in this paper. 

In the method proposed, the Green’s function in question is also obtained as the 
solution of an integral equation, the derivation of which is presented in Sec. 2. For 
simplicity’s sake, this derivation is restricted to a general two-dimensional problem 
with a single dependent variable, and deals only with the boundary value problem 
associated with Green’s function. It possesses, however, a simple physical interpretation, 
and can be readily extended to other problems. The proposed method is thus also 
applicable to three-dimensional domains, to systems of equations in several dependent 
variables, or to other boundary conditions. It can be directly applied to the determina- 
tion of Neumann’s function, and to the solution of boundary value problems in the 
theory of elasticity [3] or in the theory of heat conduction [4] (see Sec. 4). 

Two methods of solution of the basic integral equation are discussed in the paper; 
they are referred to as the series and sequence methods, respectively, according to the 
form in which the solution is most readily obtained. The series method is in effect a 
Liouville-Neumann expansion, so that the Green’s function appears in the form of an 
infinite series whose properties are well known [5]. The sequence method, on the other 
hand, has not been previously studied; it yields the Green’s function as the limit of an 
infinite sequence of functions. The convergence of this sequence is examined in this 
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paper for the case of the example of Sec. 4 (which deals with a solution of the heat- 
conduction equation), and then in Sec. 5 for the general case of Eq. (1). 

2. The basic integral equation. Consider the problem of the determination of the 
solution u(x, y) of a given linear partial differential equation for a domain D bounded 
by a closed curve C, under the boundary condition that the function u assume prescribed 
values on C. It is assumed that an integral representation of u in terms of its boundary 
values is available; for the present purposes it will be taken to be of the form: 


0G( P) a 
uQ) = [ up) SEP as, (2) 
Np 
where P and Q are points in D + C, s is the hase parameter of the boundary curve 
C, and n denotes the direction of the interior normal on C. The Green’s function G(P, Q), 
considered as a function of P with fixed Q, satisfies the conditions that 
(i) G(P, Q) = 0 for P on C 
(ii) G(P, Q) exhibits a singularity of prescribed character (depending, that is, on the 
specific differential equation in question) as P approaches Q, and is regular elsewhere 


nD+C 
(iii) G(P, Q) is a solution of the given differential equation, if considered as a function 


of P with Q fixed (and a solution of the corresponding adjoint differential equation if 
considered as a function of Q with P fixed). 

Suppose now that a function G,(P, Q) is available, which satisfies conditions (ii) 
and (iii) above, but not (i); then one may write 


G(P, Q) = Go(P, Q) + g(P, Q); (3) 


where the function g.(P, Q), considered as a function of P with Q fixed, is a regular 
solution of the given differential equation, and assumes on C known boundary values 
equal to the negative of those assumed there by G,(P, Q). Clearly then 


g(P, Q) = -|[ G(R, Q ae e.E) dsp . (3a) 
JC 


Substitution into Eq. (3) then yields the foliowing basic relation: 
aG(P, I 
AD as, (4) 





G(P, Q) = G(P, Q) -| GAR, Q) 


which may be regarded as a linear integral equation with a one kernel G, , for the 


desired Green’s function G. 
The above derivation shows that the extensions mentioned in the Introduction 


can be readily carried out; as an example, consider the boundary value problem obtained 
by prescribing the normal derivatives of u on C, rather than uw itself. It is then assumed 
that the desired solution can be written as 


AQ) = - [ N(Q, P) 2 


where the Neumann function N(P, Q) satisfies conditions (ii) and (iii) above, and the 


relation 


ou oa dsp , (5) 


ON(P, Q) 


= 2 y 5s 
np 0 for Pon c (5a) 
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in place of (i). The integral equation analogous to (4) then becomes in this case 


VP, Q) = GAP, @ + f MP, R) SHE See 2D as es 6) 


One method of solution of the above integral equations is provided by the Liouville- 
Neumann expansion. With Eq. (4) as an example, substitution of the expression for G 
given by the right-hand side of Eq. (4) into the integrand contained in that equation, 
and repetition of this operation, gives the Green’s function explicitly in the form of the 


following series: 


dGo(P, Pi) 


G(P, Q) = GP, Q) — | GAP; , Q) dsp, 
Jc Onp, 
aa IGA. .F i Gal, Pe 
+ | | GP; ,Q) — u 22) Gok P) dsp, dsp, (7) 
Jo Je ONnp, Onp, 


P F r ’ OG, dG (P2 eS 1 0G',(P3 is 0G,(P, P, 
os | | | G,(P; ,Q) —— — Py) Meet ie PY dsp, dsp, dsp, + ++> 
cwJwC dC P, Ps Ps 


Each term of this series depends only on the known function G, and on the geometry of 
the given domain, and can therefore be evaluated. The conditions for the convergence 
of expansion of this type have been extensively investigated in the past [5], and will 
therefore not be entered into here. This method of solution will be called in this paper 
the series method. 

For the integral equations arising in the present problem, it is possible to devise an 
alternative method of solution, often more advantageous, as outlined in the following 
section. 

3. Sequence method of solution. The function G, may be regarded as an initial guess 
(say the zero-th approximation) for the desired function G; other approximations may 
be obtained by considering the partial sums associated with the series of Eq. (7). The 
first approximation G, , for example, may be taken as 


GP, Q) = GAP, Q) — f Gat, ,@) SFP dy, . (8) 


This function satisfies all the conditions that were required of G, , and may therefore be 
used in its stead in Eq. (4) to obtain a new integral equation for G in terms of G, . The 
whole process may then be repeated, the next step yielding the second approximation 
G, in terms of G, , and still a new integral equation for G in terms of G, . With the notation 


1,,(P, Q) = [ GAR, Q) & oe 2A e, (8a) 
Nr 


the recurrence formula linking the (¢ + 1)th and 7th approximations is 
G,.(P, Q) = GAP, Q) — IP, Q). (8b) 
The Green’s function G and the 7th approximation are connected by the integral equation 


G(P, Q) = GAP, Q) — IP, Q), (8c) 
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where 
1(P, Q) = [ GAR, Q) ee din. (8d) 
It may be noted that the relation 
T54(P, Q) = IP, Q) — 1:P, Q) (8e) 


follows directly from Eqs. (8b) and (8c). 

Equation (8b) is the basic equation of the sequence method of solution; there still 
remains to determine, however, whether the sequence G, , G, , G, , --- converges. In the 
next section a particular example from the theory of heat conduction is solved both 
by the series and sequence methods; the methods are compared and convergence is proved 
in both cases. In Sec. 5 the case of Eq. (1) is studied, and it is found that convergence 
is assured provided the initial function G, is chosen in a prescribed manner. 

The solution procedures discussed above have been called series and sequence methods, 
depending on the form the solution most readily assumes. Clearly, however, each solution 
could be expressed in either form, and thus the two methods could be compared. If, for 





example, the third term of the sequence is referred to the original function G, , the 
following result is obtained: 
Ps IG,(P, P,) 
G.(P, Q) = GP, Q) — 2 | GAP, ,Q) EY dee, 
JC Onp, 
ref IG(P. , P:) GAP, P2) 
+2 | | GP, , Q) Moder A) GAT Fo) oe, dep, (9) 
JcoJdc Onp, Onp, 


0G,(P. ,P,) OG AP; , P2) OGA(P, Ps) 
le i — Ge Gp dip. 
Onp, Onp, Onp, 





-|[ [ [ G(P, , Q) 


e © 


and similarly for all other terms of the sequence. Comparison with the partial sums of 
the series of Eq. (7) shows that there is no direct relationship between the two sequences, 
and thus the two methods represent distinct procedures of solution. A choice of methods 
is thus available, and the advantages of one over the other may depend on the particular 
problem under consideration. It is however reasonable to expect that the sequence 
method will exhibit faster convergence since it makes use at all times of the best ap- 
proximation available rather than reverting at all times to the original guess G, . This 
expectation was fulfilled in the example of Sec. 4. where the sequence method is clearly 
preferable. The sequence method does not require the evaluation of multiple integrals; 
it furthermore provides an automatic check of the accuracy of any one approximation, 
since evaluation of the integral J;; of Eq. (8a) requires calculation of the boundary 
values of G; . 

4. Example. As an illustration of the methods of the previous sections, a simple 
problem from the one-dimensional theory of heat conduction is solved here. It is required 
to find the Green’s function for the semi-infinite solid x > 0, i.e. the temperature in 
that body at xz at a time ¢, due to a unit source released at x, at a time ¢, . This function 
will be written as G(x, x, , ¢ — t,); it is well known that its value is [4, p. 295): 


a ee even @ — 8] aa [ -et2")} 
G(x, X\ s t;) = Q[ar x(t ne tT” exp | ge t,) exp Ait 4) ; (10) 
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where «x is the thermal diffusivity; this expression will be used as a check on the present 
procedures. The temperature u(x, t) of a body bounded by a curve C, on which its value 
u(x, , t) is prescribed, is then [4, p. 294]: 


u(x,t) = « [ Ee , t) a r — )] dt; , (10a) 
70 1 Zi™zZe 





provided that the initial temperature u(x, 0) is zero. This relation takes, in the present 
problem, the place of Eq. (2); the corresponding basic integral equation analogous to 
(4) is: 
G(x. , 2, t — t) 
—1— lt, . 11 
ax, dt, ( ) 





t 
G(z, 2, ,t) = G(z,27,,) — « [ Goze » 21 » bs) 


) 


The initial guess G, will here be taken as 
Go(z,%,,t—t) = z — 7175 exp {—(z — a1)" /[4n(t — t,))}; (11a) 


that is, as the temperature in an infinite body at x at a time ¢ due to a unit source released 
at x, at time ¢, . Comparison of Eqs. (10) and (11a) shows that the difference between G 
and G, is the effect of a unit sink released at (— 2,). 

(a) Solution by sequence method. Proceeding as in Sec. 3, the first approximation 
G;, is 


2 2 


: x f' l zi x 
G(z,2,,0 =G(2z,%,) -— | ———wi eo jj -—— Chur Coo 
Na, Hs ,%1, 9) rong | [4 (¢@ — ty” exp { 4xt,  4x(t — a} i (12) 


— {ex [ ~2=20"] _ [ ~e+20"| 
” Qwxt}? \P Axl Nae Ant 


The integral in Eq. (12) may be reduced to the known form 


ew wa : } 
“172 [ exp (-- on «) dr =e’ erfe (: - x) — e” erfe (2 + x) (12a) 
a T x x 


with the aid of the transformation 





Akt’? 
= = —se 9 
f x + 4xt 7° (12b) 


Inspection of G, shows that the first approximation differs from the correct result 
by the magnitude of the sink appearing at (— 2,). Further approximations may now 
be calculated; it is found that the magnitude m; of this sink in the ith approximation 


varies as follows: 
Mm, = 0: m, = (1/2); m, = (7/8); m; = (127/128);---. (12¢) 
The general term of this sequence is 
m=1- (1/2)°°- (12d) 
and thus 


lim m,; = 1 (12e) 


io 
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and the sequence converges to the correct value. It may be noted that convergence is 
extremely rapid, and furthermore that the calculation of any one approximation requires 
only the evaluation of an integral of the same form as that of Eq. (12a). 

(b) Solution by series method. The first two terms of the series solution analogous 
to Eq. (7) are identical with those of Eq. (12); the third term has the form 


, ‘Githedatseae 7G.(t3,%2,t-h—b 
Ps [ | Go(x. , t1 , t)G(2, , 2, mf: (2 oe b=) dt, dt, . (13) 
#0 v0 X3 Xe Ie=2e 
: Leseeed 


The integrals resulting from substitution of G, into Eq. (13) reduce, with the aid of 
transformation (12b) either to the form of Eq. (12a) or to the form 





a= a. 


4 ‘ » aie 
=o | 7 exp (=? _ ) dr = (a+ le”. (13a) 
7 J0 \ T 


The final result is 


1 f (2 
G(x, , ) = sia exp | "| (13b) 
3b) 


2% x + 2,)* || 
_ (0 Tet 4 + ++) exp | —@ +20" ]} 


which shows that the strength of the sink at (— 2,) increases geometrically with each 
term. Convergence is thus assured, though it is not as rapid as that of Eq. (12c). 

5. Convergence of the sequence method in the special case of Eq. (1). Equation (1) 
was chosen for this detailed investigation of the sequence method because of its im- 
portance in theoretical physics, and because its properties have been widely studied 
(see, for example, the comprehensive treatise of Bergman and Schiffer [1], which will 
be taken as a general reference for this section). 

The singularity required of the Green’s function G(P, Q) of Eq. (1) as P approaches 
Q is a logarithmic infinity; it will further be specified that 


. G(P 
lim | aG(P, Q) dsp = —1, (14) 


Onp 


r—0 


where c, is a circle of radius r around Q. Functions which satisfy Eq. (14) and conditions 
(ii) and (iii) of Sec. 2 will be called fundamental singularities; between any two such 
functions, say S,(P, Q) and S,.(P, Q) the relation 

; ace ' aS.(R, P - aS(R, 

S\(P, Q) = S.(Q, P) = f | suc, Q) Sa, P) ose SAR, FD 988, dsp (14a) 

Cc ONnR Onr 

must hold [1, p. 270]. Green’s function is of course an example of a fundamental singu- 
larity, and it follows from Eq. (14a) by setting S,(P, Q) = G(P, Q) and S,(Q, P) = 
G(Q, P) that it is symmetric, in the case of the present (self-adjoint) differential equation, 
i.e. 


G(P, Q) = GQ, P). (14b) 


It is interesting to note that in the case of Eq. (1) the basic integral equation of the 
present method, Eq. (4), follows directly from Eqs. (14a) and (14b) by setting S,(P, Q) = 
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G,(P, Q) and S.(Q, P) = G(Q, P) and transposing terms. Equation (6) can be similarly 
derived. 

In order to examine the convergence of the sequence G, , G, , G, , --- defined by 
Eq. (8b), the following inductive reasoning will be employed. It will be first assumed 
that a particular member of the sequence, say G, , has the following properties (in addition 
to those previously prescribed for all members of the sequence): 


(i) G, is continuous and twice continuously differentiable in a domain D*, bounded 
by a curve C*, which contains the given domain D; except possibly on the bounding 
curve C of D, 

(ii) G,(P, Q) = 0 for P on C%*, (15) 

(iii) G.(P, Q) = G,(Q, P). (15a) 


It will then be shown (a) that the function G,,, also enjoys properties (i), (ii) and (iii) 
above, and (b) that under these conditions the sequence G, , Gi, , Gis2 , +++ converges. 
If then the initial function G, is chosen so as to satisfy the above conditions, convergence 
of the sequence G, , G, , G, , --- will clearly be assured. 

To prove statement (a) first, consider the identity 


IP, -[ GAR, Q) GRD dsp 
ONE 


(16) 


ee P) 


= Ep» (G(R, P); G(R, Q} + ‘ G,(R, Q) dts 


which follows from Green’s theorem [1, p. 260]. The last integral of Eq. (16) vanishes 
because of Eq. (15), and the energy integral £, extended over the domain indicated in 
the subscript, is defined as follows: 


E{u(R, B);(R, Q)} = [f (ae? av(R, Q) 


OXp OXR (16a) 
u(R, P) WR, QD) 1 opyu(R, Pr(R, @} des de. 
OYR YR 


It is thus clear that J,, is symmetrical, i.e. 
Tu(P, Q) = In(Q, P) (17a) 
and that therefore 


aG.(R, Q) 


AGAR, ( : 
— D ase = | GAP, p) ee dee (17) 


IP, Q) -[ G(R, P) 
and finally that 
Tn(P,Q) =0 for PonC*. (17¢) 


Inspection of Eqs. (17a, b, c) proves statement (a) above. 
Before proceeding with the proof of statement (b), it is convenient to establish the 
following inequality 


> In(P; , Pix; 2 0 (18) 


t-1 j=1 
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for any choice of points P; and P; in D and of real numbers 2; and zx; . The quadratic 
form in question is in fact identical with the expression 


N N 
zz: bo E'p--p{G.(R, P;:); G(R, P;)} 2,2, 
5 7-1 (18a) 


N N 
= Emo > GAR, P;)z; } 2 GAR, P)a)} Fad 0. 


It is furthermore shown in [1, p. 279] that if the function q(z, y) of Eq. (1) is analytic 
in x and y throughout D + C, then the equality of Eq. (18) can hold (the trivial case 
x, = 0 excepted) only if 

I,, = 0. (18b) 


Consider now the quantity J;,, ; it may be written [see Eqs. (8c-e)] as: 
, ; P 0 ; 
Ii(P, Q) = 1(P,Q — LP, = | GAR, P) = (G(R, ® — GAR, Q)] se 
J , Onp (19) 


al 1, Y) , 
a 


lI 


-[ 18,P) sp = Ep{IA(R, P); 1(R, Q)}. 
From the latter energy integral, by means of a reasoning analogous to that which led 
from Eq. (18a) to inequality (18), it follows that the quadratic form 


N N 
> SLAP: , Pra; k¥0 


#=1 j=l 


is positive definite. It follows from Eqs. (8e) and (18) that the quadratic form associated 
with (J, — J,.,) is also positive definite, and that therefore the relation 


> y LP. , P;)x;2; > > fs (P;,P)22;>0; kx#0 (19a) 


i=l 
holds. This implies that the sequence of the quadratic forms associated with J, is a 
seat decreasing sequence and therefore converges (and does so, in fact, mono- 
tonically). According to Cauchy’s convergence theorem, then, corresponding to any 
given positive number e, however small, there exists a number m(m > k) such that the 
inequality 


N N V N 
b >. Li, Para, | ~_ }> p ag OE A Pa; | <e€ (20) 
#=1 j=1 i=l j=l 


holds for all positive integral values of p. With the aid of Eq. (8c) this may be written as 


N 
i=1 


N 


[7.AP; » Pi) — InsP; , P;) 2:2; 
=1 


? 


(20a) 
N V 

> > 1G,(P:; ,P)) — Gui AP; , Pia; < €. 

i=l j=l 


This inequality must hold for arbitrarily selected x,’s; if they are all chosen zero except 
x, = 1, say, it follows that 


G.(P, , P:) — G.sAP, ,P) <<, (21a) 
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where P,; may designate any one point in the domain considered. Similarly, one may 
write 
Gn(P2 ,P2) — Gn+o(P2 , P2) < € (21b) 


for another point P, . The quadratic forms in (20a) are however positive definite, so 
that their principal minors of order 2 must be positive; in other words: 


| Gu(P: , Ps) — Gnsp(P: , Ps) | ina 
< (G,.(P fs — Gmi(Pi » P,))]'1G,(P2 , Pe) — Gas(P2 PI” < €. 


Inequalities (21a) and (21c) then show that the conditions of Cauchy’s theorem are 
everywhere satisfied by the sequence G, , G,,, , «++ , and that therefore this sequence 
converges. 

There now remains to select the initial function G, so as to satisfy the conditions 
listed earlier in this section. It is clear that this will be the case if G, is chosen as the 
Green’s function of a domain D* which includes D. This is probably the simplest choice, 
though by no means the only one; for example, the Neumann function of D* will also 
be a suitable selection, though probably less efficient in view of the relation [1, p. 383]: 


N(P, Q) = GP, Q) (22) 


for any domain. The proof in this case is analogous to that just presented. 

The above proof shows that convergence is assured, but it indicates that the limit 
reached need not necessarily be, from a general standpoint, the desired Green’s func- 
tion; for example, Neumann’s function also satisfies Eq. (18b). This situation arises 
from the fact that in the set up of the sequence method, only the first two terms of 
the series expansion of Eq. (7) were employed; on the other hand, it appears impractical 
to apply the method on the basis of additional terms. This drawback of the sequence 
method is not serious, however, from the viewpoint of a calculation procedure. It is 
usually unlikely that any function but Green’s will result, since the portion of the ex- 
pansion (7) used pertains to Green’s and no other function; and in the example of the 
previous section no such difficulty was encountered. It was mentioned in Sec. 3, moreover, 
that at each step of the sequence method a running check is automatically kept on the 
improvement obtained. Thus, if it should happen that at any one stage of the calcula- 
tions convergence appears unsatisfactory, the work already performed is not wasted; 
the integral equation may be set up with the last (and still best) approximation available, 
and the series method applied to that. 
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BOOK REVIEWS 


Electrodynamics. Lectures on Theoretical Physics, Vol. III. By Arnold Sommerfeld. 
Translated by Edward G. Ramberg. Academic Press, Inc., New York, 1952. xiii + 
371 pp. $6.80. 


Following a well written historical review, the author introduces the integral form of Maxwell’s 
equations on an axiomatic basis and obtains the transition equations at a boundary between two different 
media. After Maxwell’s equations have been put into differential form the effect of the properties of 
the media are considered and with the concepts of energy flow the author introduces Poynting’s vector. 
While the rationalized MKSQ system is used there are long discussions on the choice and number of 
fundamental units. This is supplemented by an excellent table at the end of the book giving the dimen- 
sions of symbols used and numerical values of certain constants in a logical presentation, i.e., c = 
velocity of light, (measurement); 4) = permeability of vacuum, (definition); ¢g = dielectric const. 
of vacuum, (consequence). 

In Part II, which comprises half of the book, the usual complement of subjects of classical electro- 
dynamics is treated with the welcome addition of remarks on ferromagnetism. Here the author comments 
on the Weiss domains, electron spin as an elementary magnet and emphasizes that “‘. . . ferromagnetism 
is not based on Maxwell’s phenomenological theory, but on the more profound laws of atomic physics 
and on the statistical behavior of the electron.” 

In Part III the theories of relativity and electrons are covered and numerous references are given 
to the original papers on these subjects. The author clearly points out the limitations of the electro- 
dynamic theory of the electron and restates that, at present, quantum theory is the only way in which 
a definite electron spin and magnetic moment can be defined for the electron. 

The last portion of the book, Part IV, considers Maxwell’s theory for moving bodies and brief 
comments are made concerning the general theory of relativity on the unified theory of gravitation. 

This volume, a valuable reference for those more experienced with the material, would be well- 
suited for use as a text in an advanced course in electrodynamics. It does have the disadvantage that 
all the problems are worked out, requiring an additional source of problem material. 

SHevpon L. Levy 


Advanced statistical methods in biometric research. By C. Radhakrishna Rao. John 
Wiley & Sons, Inc., New York and Chapman & Hall, Ltd., London, 1952. xvii + 
390 pp. $7.50. 


This is a thorough account of the methods and the theory of the analysis of variance, multivariate 
analysis and related topics. The first two chapters collect the necessary mathematical background. 
The first, the elements of matrix algebra and the second, distribution theory with particular emphasis 
on the normal distribution and the distributions arising in the analysis of variance and in multivariate 
analysis. The next chapter deals with linear estimation and tests of linear hypothesis. This is followed 
by a discussion of estimation in general and in particular the theory of the maximum likelihood estimate. 
Fisher’s scoring method on iterative procedure for solving the maximum likelihood equations in more 
complicated cases is explained and exemplified but the question when this procedure converges is not 
discussed. After a discussion of large sample theory there follow three chapters on multivariate analysis, 
starting with tests for homogeneity of variances and correlation coefficients. After a discussion of dis- 
criminant functions and of Wilks’ criterion the author devotes a chapter to the application of discriminate 
functions to classificatory problems. The last chapter discusses measurements of racial difference, in 
particular Mahalanobis’ generalized distance and Karl Pearson’s coefficient of racial likeness. 

The applications and examples are drawn from biology and related fields but the book does not 
overemphasize these applications and will be of equal value in every field of research, where the need 
for the use of statistical methods arises. 


H. B. Mann 


(Continued on p. 266) 
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ELASTIC-PLASTIC TORSION OF CIRCULAR RING SECTORS* 


BY 
W. FREIBERGER 


Brown University 


Abstract. The solution is presented to the problem of uniform torsion of circular 
ring sectors with circular cross section under the assumption of perfect plasticity. The 
elastic-plastic problem is solved by a semi-inverse method. When the entire ring is 
plastic a discontinuity of stress appears which may be regarded as the limiting case of 
an elastic core; such discontinuities have recently been discussed in the literature [1]. 
The solution for the discontinuous fully plastic stress distribution is exact, that for the 
elastic-plastic case approximate in the sense that it is found exactly for cross sections 
differing slightly from the circular. This difference is negligible for ratios of ring radius 
R to cross section radius p occurring in practical applications to helical springs of small 
pitch. 

1. The completely plastic problem. The stress system being independent of the 
meridional coordinate @ the equations of equilibrium in cylindrical polar coordinates 
reduce to the single equation 





a(r’r,) + a(r*7,) 7 


or dz 0, (1) 


where r, and 7, are the only non-vanishing stresses, viz. the shearing stresses transmitted 
across meridional planes. The yield conditions of Tresca or von Mises require that 
2 2 2 ‘ 
Ty + T; _ k ? (2) 
where k is the yield stress in simple shear assumed to be a constant throughout the body. 
The problem is completed by the boundary condition that the resultant shearing stress 


with the components 7, and r, be tangential to the contour of the cross section. 
To solve Eqs. (1) and (2) set 


t, = k cosy, tT, = ksin y, (3) 
where y is the angle between the stress vector and the positive r-direction; hence 

Oy Oy 2 me 

ar cot y ae : cot y = 0. (4) 


This equation can be solved by the method of characteristics. The characteristic equa- 
tions are 
dr ; dz dy 2 
= sin — = —cos y, —- = — cos yp. 5 
ds v ds v ds r v (5) 
It follows from (3) and (5) that at any point the resultant shearing stress vector is 
normal to the characteristic through that point. In particular, the characteristics are 
normal to the stress-free boundary. 
Integration of the first and second of Eqs. (5) gives 
r’ cos y =, (6) 


*Received Sept. 21, 1955. 
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where c is a constant along a given characteristic. From (6) and the first two Eqs. (5) 
it follows that 


dz { (on Jnl /2\4 -1/2 
pn eae = Fign ) =— Ty : 
or 
z = F(c/2)'"F(45°, arcosc’’?/r) + constant, (7) 


where F(a, ¢) is the incomplete elliptic integral of the first kind, in the notation of 
Ref. [2]. The constant in Eq. (7) is determined by the boundary condition along the 
stress-free surface; if 6 is the angle between the positive r-direction and the radius 
vector from the center of the cross section to a boundary point P, we have 7, cos 8B + 
7, sin 8 = 0, or, by Eq. (3), cos (W — 8) = 0 and hence 


sin ¥y = cos at P. (8) 
Moreover, at P, 
z= psinB, r=R-+ peosB. (9) 


Hence, from Eqs. (6), (8), (9) the constant c in Eq. (7) has the value 

c = (R + p cos B)’ sin B (10) 
for the characteristic through boundary point P(p sin 8, R + p cos 8). With this value 
the equation of the characteristic finally becomes 


; and eee een (R cos 8) sin’ 
z/p = sin B + 2” ( + cos 8} sin'”’ 6 F(45°, arcos emia © 
\p \ r 





(11) 
— F(45°, arcos sin’”’ B) 


The value of ¥, and hence that of the shearing stress components, on any point along 
the characteristic (11) is found from the equation 


R cos B\* . ; 
cos Y = (Rte cose) sin p. (12) 


r 


R_ 
po! 
choracteristics 








direction of 
shegring stress 
vector ® 














Fic. 1. Fully plastic stress distribution: characteristics and line of discontinuity for two ring sizes. 


The characteristics intersect on a line of stress discontinuity. Equilibrium allows only 
the tangential component of the shearing stress to be discontinuous. The line of stress 


discontinuity must therefore bisect the angle formed by the characteristics meeting on 
= 4. 


this line. Figure 1 shows this situation for the rings specified by R/p = 1 and R/p 
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It can be shown (see [3], [4] and [5]) that the resultant of the considered stress system 
is indeed a force acting along the axis of the ring; the expression for this resultant is 


Z 7 R * sin? B 1/2 (2) (2) 
kp ll {1 i (: la a) sth . p . p/’ 


with the integration carried out over the cross section. 
2. Velocities. The velocity-strains for the present problem obey the relations 





Ov, Ov, OV, Ov, Ov 
= 0, —_=-(0 : at ) = , cml —- 
or Oz 00 + 0 or Oz 


l (. dv, Os ) a 1 (2 1 *.) os 
2\r 20.7 or r = Pe 2 ie 00 = Xr; 


which are satisfied by velocities of the form 


v, = 0, = Ke, . = g(r, 2), (13) 


= 0 


where K is a constant and ¢(r, z) denotes the warping function. Thus, 
dp _ ¢ 


re) 
— —£ Dr, dg | K 


— 
or r dz r mts (14) 


From (3) and (14) the differential equation of the warping function is seen to be 


d(y/K) _ cot ¥ + o/K 


O(e ‘K) = 
4 0z r 


Ls cot y 


rs = 0. (15) 


Equation (15) has the same characteristics in the r, z plane as Eq. (4), and the variation 
of ¢ along the characteristics is given by 
die/K) _ o/K _ _lde 


dr r r dr 





By integration one obtains: 


4 _ r/p - 1/2 1/2 { r R/p + cosB . i/2 
+ sin B42 B( 45° arcos ———————— sin 
K "R p + cos ps af : r/p B 


ye / 
i 2-VF( 45°, arcos sie + ose sin” a)} a 


where F and E are the incomplete elliptic integrals of first and second kind, respectively, 
and the constant C is determined from the boundary condition on the line of discon- 
tinuity. Since the line of discontinuity is the limit of a narrow elastic core, the shear 
rate of the assumed rigid, plastic material must vanish along the line of discontinuity. 
Thus, 


dg ) dg K 
or r Oz r 


and hence 
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Since the circumferential velocity component v» is only determined to within an arbitrary 
additive constant, which corresponds to a rigid body rotation of the ring about the z 
axis, the function “” may be given the value zero at some point on the discontinuity. 
Numerical integration of (15) then furnishes ¢ along the discontinuity and so the re- 
quired initial condition for integration along the characteristics originating on the 
discontinuity. Examples for this procedure are provided in Refs. [4] and [5]. 

3. The elastic-plastic problem. The solution to the elastic-plastic problem of a 
ring with approximately circular cross section will be obtained by an inverse method 
first proposed by Sokolowsky for the torsion of an elliptic cylinder [6]. 

It can be shown (see [7]) that the elastic torsion problem of circular rings requires 


te ee (17) 


\ , oon r 


the solution of the equation 


Here, ® is a stress function in terms of which the stress components are expressed by 


A 0® A 0® ’ 
t, = T= -3 , (18) 


r 02 ‘ Yr or 
where A is a constant affecting the absolute values but not the distribution of the stresses. 
For a given ring, this constant is proportional to the applied force. 
Equation (17) can be written 
edb 306 , dS 
ay ae Sa py ae ae 


= —z l. (19) 
or fr or Oz ‘ 


If a solution of Eq. (19) can be found giving a stress distribution such that the yield 
condition Eq. (2) is satisfied along a more or less circular curve, then this curve may be 
taken as the elastic-plastic interface and used as the starting curve for the construction 
of the characteristics emanating from it. This is done by evaluating the ratio r,/r, along 
the curve from the elastic solution and noting that the characteristics will be perpendicular 
to the direction of the vector (r, , 7,). Any orthogonal trajectory of the characteristics 
can be taken as the boundary of a cross section, for which an elastic, plastic stress dis- 
tribution has thus been obtained. If two or more solutions of Eq. (17) give practically 
the same cross section, the progress of the plastic region with increasing external force 


can be traced. The stress function 
® = (r? — a’)? + D2’ (20) 
is found to satisfy Eq. (19) if the constants a and b are related by the equation 


— 8a’ 
b= 2 i (21) 


9 


Substituting Eq. (20) in Eq. (2) it follows that possible elastic-plastic boundaries are 
given by 


cz = ra — (r’ — B)’), (22, 


where 


2 2 
c’? = b*/4, a‘ = i a + sa’) / 644°, B= (*, + sa’) /'s, (23) 
16A 4A 
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and a and b are related by Eq. (21). Clearly, 8* — a* = a‘. Figure 2 shows two elastic- 
plastic interfaces; they are found from Eqs. (22) and (23) by giving the parameters the 
values a = 1.5, 8 = 3 (interface 1) and a = 1.1, 8 = 3 (interface 2), respectively. By 
drawing the characteristics based on these curves it is seen that the cross sectional 
boundary R/p = 6 in Fig. 2 is approximately an orthogonal trajectory to the families 


based on either. 












Zz 
zf 
characteristics based on interface 2 
Sh , 
cross-section boundary 
4 of ri R = 
ng Pp 6 
.3f elastic-plastic interface | 
ol .. elastic-plastic interface 2 
| b 
eae nicer adh Sas sae al otin , 3 
5.0 34 58 62 66 “>? 


Fic. 2. Elastic-plastic solution. Showing boundary of ring considered and two elastic-plastic interfaces. 


It remains to evaluate the resultant force corresponding to the fully elastic ring, the 
fully plastic ring, and the ring for which the elastic-plastic boundary has advanced to 
either of the two positions of Fig. 2. If the resultant force is denoted by Z, 

Z 1 ff 


4-1 | r, dr/p) de/p) 
cp C Je 


r | fle 
- /| sin ¥ d(r/p) d(z/p) — i” I a a d(r/p) d(z/p) (24) 


il I/ sin wy d(r/p) d(z/p) — 4 A I (" -- ) d(r/p) d(z/p), 


where the first integral is taken over the plastic, the second over the elastic part of the 
cross section. Both integrals are best evaluated numerically by dividing the cross section 
into a net of squares and summing the values of the integrands at the centers of these 


squares. 
Fora = 15,8 = 3 one finds a’ = 8.71, A/k = .329, 
fora = 1.1,8 = 3 a” = 8.96, A/k = .625. 
The results of the integrations are the following: 
ting fully elastic, with yield value just reached 
at point of maximum stress: Z/kp? = .2 
Ring with elastic-plastic boundary 1: Z/kp’ = .22 (25) 
Ring with elastic-plastic boundary 2: Z/kp’ = .305 
Ring fully plastic: Z/kp’ = .328 












264 W. FREIBERGER [Vol. XIV, No. 3 


It is shown in Ref. [7] that, approximately, 
Z/kp” = (p/R)*(uKx/2p’) (26) 


for the elastic ring; here, K is the constant in Eq. (13). The first three equations of Eq. 
(25) substituted into Eq. (26) with the value of p/R appropriate to the size of the elastic 
core in question give three points on the graph (external force) vs. (axial velocity), and 
the fourth equation of Eq. (25) gives the value of the asymptote for K — ~. This graph 
is plotted in Fig. 3. 











! i i 1 i iL i i =. 
100 200 400 600 800 1000 1200 2k 


Fic. 3. Resultant force vs. axial velocity. Resultant force: 7 


Axial velocity: vz = K@. 
A: End of purely elastic deformation 
B: Elastic-plastic interface 1 of Fig. 4 
C: Elastic-plastic interface 2 of Fig. 4 
D: Resultant force for completely plastic ring 
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Simultaneous linear equations and the determination of eigenvalues. Edited by L. J. Paige 
and Olga Taussky. Issued by the National Bureau of Standards, Applied Mathe- 
matics Series 29, 1953. iv + 126 pp. $1.50. 

This volume contains the majority of the reports presented at a Symposium which was held August 
23-25, 1951, in Los Angeles, under the sponsorship of the National Bureau of Standards, in cooperation 
with the Office of Naval Research. The principal topic is that of approximate and numerical methods of 
solution and of the determination of eigenvalues. The volume contains the following papers: Tentative 
classification of methods and bibliography on solving systems of linear equations, by George E. Forsythe; 
Simultaneous systems of equations, by A. M. Ostrowski; The geometry of some iterative methods of 
solving linear systems, by Alston S. Householder; Solutions of linear systems of equations on a relay 
machine, by Carl-Erik Fréberg; Some special methods of relaxation technique, by Eduard Stiefel; 
Errors of matrix computations, by Paul S. Dwyer; Rapidly converging iterative methods for solving 
linear equations, by J. Barkley Rosser; Some problems in aerodynamics and structural engineering 
related to eigenvalues, by R. A. Frazer; Inclusion theorems for eigenvalues, by H. Wielandt; On general 
computation methods for eigenvalues and eigenfunctions, by Gaetano Fichera; Variational] methods 
for the approximation and exact computation of eigenvalues, by Alexander Weinstein; Determination 
of eigenvalues and eigenvectors of matrices, by Magnus R. Hestenes; New results in the perturbation 
theory of eigenvalue problems, by F. Rellich; Bounds for characteristic roots of matrices, by Alfred 
Brauer; Matrix inversion and solution of simultaneous linear algebraic equations with the IBM 604 
electronic calculating punch, by George W. Petrie, III; Experiments on the inversion of a 16 X 16 matrix, 
by John Todd; A method of computing eigenvalues and eigenvectors suggested by classical results on 
symmetric matrices, by Wallace Givens; Computations relating to inverse matrices, by Jack Sherman; 
Results of recent experiments in the analysis of periods carried out in the Istituto Nazionale per le Appli- 


cazioni del Calcolo, by Gaetano Fichera. 
G. W. Morcan 


Handbook of probability and statistics with tables. By Richard 8. Burington and Donald 
C. May, Jr. Handbook Publishers, Inc., Sandusky, Ohio, 1953. ix + 332 pp. $4.50. 
“It is the purpose of this book to provide a small convenient size pocket handbook of probability 

and statistics sufficiently comprehensive to fill a wide variety of needs, ...”. This statement of the 

authors seems to the reviewer to be a modest statement of the book’s accomplishments. The summaries 
of frequency and probability distributions, generating and characteristic functions, regression theory 
and time series, sampling distributions, statistical inference, analysis of variance, finite differences and 
interpolation, sequential analysis and quality control constitute the bulk of this handbook. In addition 
there are various tables of special distributions as well as a rather complete reference list. While this 


book claims only to be a handbook it could be used as a text as well. 
Roun TRUELL 


Partial differential equations in engineering problems. By Kenneth 8S. Miller. Prentice- 

Hall, Inc., New York, 1953. viii + 254 pp. $4.75. 

This text presents in a simple and straightforward manner the classical techniques of Fourier series, 
integrals, and transforms and applies them to find solutions of boundary value problems associated with 
the common partial differentia] equations of mathematical physics. It is written on the level of a student 
who has an elementary knowledge of partial differentiation and ordinary differential equations. Every 
significant concept and technique is illustrated by a specific problem worked out in complete detail 


and each chapter ends with an extensive set of problems. 
P. CHIARULLI 
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CHARACTERISTIC VALUES OF ARBITRARY MATRICES* 


BY 
MARK LOTKIN 


Avco Advanced Development Division, Stratford, Conn. 


1. Introduction. The determination of the characteristic values of matrices is of 
paramount importance in many branches of the applied sciences, and consequently, 
much attention is being given to the improvement of existing relevant methods as well 
as to the development of new ones. 

The method described below for the determination of the characteristic values of 
arbitrary matrices possessing complex elements, aside from being new, is believed to be 
of promise. It is based on a theorem of I. Schur which asserts that any matrix of complex 
elements can be reduced to triangular form by means of sequence of unitary trans- 
formations which generate the characteristic values of the matrix along the diagonal 
of the triangularized form; see, for example, MacDuffee [1]. In applying this theorem, 
as described in detail below, the triangular part of the matrix initially possessing the 
lower norm is chosen for annihilation. Then an arbitrary element in that triangle, e.g. 
an element of above-average absolute value, is selected as the pivot of a certain unitary 
transformation which is determined in such a manner as to reduce the total norm of 
the triangle. The construction of the unitary transforms having this norm reducing 
property is founded on the real roots of associated cubic polynomials with properly 
determined coefficients. 

The entire process is iterated until the norm has been decreased to a desired tolerance. 

In the following sections, then, the effect of the transformation upon the matrix is 
investigated, and the transformed elements are exhibited individually. In Sec. 3 the 
behavior of hermitian and skew-hermitian matrices is considered briefly. The con- 
struction of the norm reducing cubic is discussed in Sec. 4. The remaining two sections 
are devoted to the applicability of the proposed method to high speed computing 
machinery, and to a simple example illustrating the procedure. 

2. The unitary transformations. Let us consider the n-dimensional matrix of com- 


plex elements 





*Received October 20, 1955. 
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Here the elements a, b, c, d, are, respectively, in the (ii), (ij), (ji), and (jj) position of A. 
It is assumed here that the “upper” triangle of A is to be annihilated. The element b 
will be called the “‘pivot’”’ of the transformation. As shown below the matrix A can be 
reduced to triangular form by an infinite sequence of unitary transformations 7, , 
p = 0, 1, 2, --- , so that 
Aug @ TEATS , A, =A, 

where 7% = T/ , denotes the conjugate transpose of 7, . As p grows beyond bound the 
diagonal elements of A, then should approach the characteristic values of A. We choose 
the transforms to be of the type 


| ‘ 0 
t —(l-—r)” - |] 
| 
T = l a (2 
(1 —r’)'” t 
_0 . ] 
with 
t = r exp (76), r real. 
Thus T contains ¢, — (1 — r*)'””, (1 — r’)’”, and the complex conjugate 7 of ¢ in the posi- 
tions corresponding to a, b, c, d, and agrees with the unit matrix everywhere else.’ The 
number ¢ will be appropriately determined later. Clearly 7 is unitary: T~' = 7*, and 
a , 0 | 
i (l—r)'”” 
| 
T* | (3) 
—(1-—?r)” t 
0 F ‘ ; 1 | 





1This transformation was suggested by S. Pines, Republic Aviation Corporation. 
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Let us now form 
A, = TAT = (a{}). 
For the sake of simplicity the subscripts ‘0’”’ have been deleted here. In calculating A, 
it is convenient to partition T as 
| Ores.ieian = Og-1 0-3 
PT Bicccan. ~Gesis “CGeistis (4) 
\ tes | (Os ects | 


with J, denoting a unit matrix of k dimensions, 0,, denoting a rectangular matrix of 
zeros having r rows and s columns, and, finally, S denoting a square matrix of 7 — i + 1 
dimensions. Similarly, let us denote the corresponding parts of A by 


(B C D 
A= le F Gi, (5) 
| H JK | 
having, respectively, the dimensions 
B: @—1)xX(@—1) C: @—1)xX(j—ti+1) D: @—1)X(m-—j)) 


E:(j-it+)xXG@—-) F:G-it+)XG-it+d) G:G-it)xa-jJjy 
H: (n—j)X@—1) J: (n-j)X(j—-t+) K: (n — j) X (n— 9). 


It is then observed that 

B Cs D | 

A,=|S"'E S'FS S"G |. (6) 
H JS K | 

The elements of the central part S~' F'S of A, are thus seen to be constituted as follows: 
Let 


R=(1-9r)'-r", z = R exp (i8). (7) 
Then 
a, = r'[a + R’d + zc + 2b], (8) 
b, = F[b — 2c + 2(d — a)], (9) 
ce, = tle — 2b + Xd — a)], (10) 
d, = r'[d + R’a — zc — 2b], (11) 


for the corner elements of S”'FS. 
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The other transformed elements occur only in the 7th and jth rows and columns of the 
affected components of A, , and are obtained as follows: 


(i) ith row: an = U ua t+ 2a,), (12) 
(ii) jth row: a}, = taj, — 2a;), (13) 
(iii) 7th column: ay; = tay; + 2a,,;), (14) 
(iv) jth column: ay; = Ua; — 2a,;). (15) 
The range of k is given by: 

For CS: Formulas (14) and (15) forJ,:1<k<i-—1, 

For JS: Formulas (14) and (15) for/J,:j7 +1<k <n, 

For SE: Formulas (12) and (13) for 7, , 

For SG Formulas (12) and (13) for J, , 





For S-'FS: Formulas (12) through (15) for 7; :7 +1<k<j-—1. 


In addition to the relationships derived above it.must also be remembered that unitary 
transformations leave trace and norm invariant: 


>> | aS? |? = constant. (16) 


2 
t,j3=1 


>> a {? = constant, 
1 


The determination of suitable choices of real-valued R, @ is discussed in Sec. 4. Once 
R, 6 are known the quantities 7, ¢ are obtained as 
r= +(1 +R’), t = r exp (76). (17) 

The sign of r does not affect the new corner elements a, , b; , c, , d; . However, the signs 
of the other elements of A, are changed if the sign of r is changed. To fix matters it is 
advisable to assign to r the same sign as is possessed by R. It is seen that |r| < 1 for 
R ~ 0, while for R = 0 also z = 0, and r = 1. In the latter case the transformation 
leaves the diagonal elements unchanged, rendering zero values of FR of no interest. 

3. Hermitian and skew hermitian matrices. Let us pause to consider the effect of 
the transformation 7’ upon hermitian and skew hermitian matrices. 

A. The matrix A is hermitian. Here a,; = 4;, for all 7 and k. Thus the diagonal elements 
are real, and c = 6. Equation (8) becomes 


a, = r'[a + R’d + zb + 2b], 
so that a, is real if R is real, which will be assumed in the following. Equation (11), or 


the fact that a, + d, = a + d, shows the diagonal element d, also to be real. Equations 
(9) and (10) reveal that c, = 6, , and the other transformations (12) through (15) indicate 


(1 =(i) a] . . . °,° 
that a,;\ = 4d;, , in general. Thus A, , and, similarly, A, is also hermitian. 
B. The matrix A is skew hermitian. In this case a,; = — 4, fort, k = 1, 2, ---,n, 
so that a;; = 0 for all 7, c = —b. Equation (8) leads to 
. = 
a, = r lec — zc], 
making a, a pure imaginary number, or zero. Since a, + d, = 0,d, = — a, is also purely 
imaginary. Further, it turns out that c, = — 6, , and aj;’ = — dai,’, i ¥ k, in general. 
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Continuing with the second iteration we find that 
a, = r'[(1 — R’)a, + 2b, — 2b,], 


d, = —@Q.. 


~(2) 


Thus a, , d, are also pure imaginaries. Next, c, = — 6, , and a;;’ = — di,’ . The successive 
transformed matrices A, all show the same behavior, namely, that the skew hermitian 
character is preserved except for the elements of the diagonal which become pure 
imaginaries. 

The preservation of the hermitian and quasi skew hermitian nature of the trans- 
formed matrices is clearly of importance in that it essentially halves the amount of 
computational labor involved in the entire calculation. 

As has been shown by the author elsewhere, hermitian and skew hermitian matrices 
can be diagonalized by unitary transformations of the type (2) by the annihilation of 
the pivot. 

4. Construction of the norm reducing cubic. It is desired to construct the variable ¢ 
in the transformation matrix T in such a manner that the sum of the squares of the 
absolute values of all the super-diagonal elements in the transformed matrix A is de- 
creased. Now this sum, for the affected elements of the original matrix A, is 


M = Di(jau + lay?) + D (au |? + | ae |) 
1 2 
+ > (| ax > + | a; ") + |b)? (18) 
3 
=> + D+ D+/5)', 
1 2 3 
and the summations are to be extended over J; for > ae , 7 = 1, 2, 3, with the intervals 


I; as defined in Sec. 2. For the transformed matrix A, the corresponding sum is 


(1) 


M= D+ 4+ L410, (19) 


where the superscripts ‘‘1’’ indicate that the a,, are to be replaced by the corresponding 
_. 
Now the transformation formulas, together with the fact that r? (1 + R’) = 1, lead to 


(1 


Dar V+rR’ DV (lan |? + | a: |’) 
3 3 3 
Ho? De [leajdin + 2djQin) — (20,44; + 2diay;)). 


Consequently, 
M,-—- M=(r? —1) + 2rR’ YD (| an |? + | ax: |) 
+ 2°R > [| aj | | an | cos (6 + an — an) 
- | Aes | | Ay; | cos (0 + ax; — onj)] + | by ? - |b |’, 


where in general a,; = | a;; | exp (ta,;). 
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Next it is found that 


lb, P=r*{{bP?+R*tlce|? +R? |d—al|’? — 2R’|b| |c| cos (26 — B+ y) 


+ 2R|b| [| d|cos(@+ 6 — B) — |a| cos(@ +a — B)] (20) 
— 2R° |c|[|d|cos(6+ y — 6) — |a| cos(0@+ 7 —a)}}, 
where again a = | a| exp (za), ete. 


For |7| < 1, then, 
M,-—M < RF(R, 9), 


where 
F(R, 0) = C.R°+C.R°+CR+C, (21) 
C= |e (22) 
C, = —2!/c|[\d|ecos(@+y— 6) —/|a|cos(@+ 7 — a)] (23) 
C,=!d-—a|\’—2|b c|cos(20—B+y) + > a a: |° + | aj |°) (24) 
C, = 21/6! [\d!|ecos(@+ 6 — B) —!a|ecos(@+a — £)] (25) 
+2 Pd a, a;, | cos (0 + aj, — ayy) — | a; a,; | cos (6 + a,; — a,;)). 


For any value of @ the cubic polynomial F(R, @) has at least one real root; the real 


root of largest absolute value leads to an 
r = (R’+ 1)” 


of least value such that M, < M. Thus the system 


F(R, 0) = 0, OF(R, 0)/00 = 0 (26) 
will produce the z = R exp (76) required in the transformation. 

It is sufficient to restrict the values of 6 to the interval (0, x). If, namely, F is a real 
root of F(R, 6) = 0 associated with a value of @ e (0, 27), then, if we put 6, = 6 — 7, 
R, = — R, because of 

C,(6,) = —C,(8), C,(0,) = —C;(8), 
it follows that F(R, , 6,;) = 0. However, 
z, = R, exp (10,) = R exp (70) = z, 
and, if sgn r = sgn R, also 


t; = —r exp (70,) = r exp (20) = l. 


The value of z as determined by the system (26) will produce the greatest possible 
reduction per iteration in the superdiagonal norm. Any other number z obtained by 
solving F(R, 0) = 0 alone for arbitrary @ e ‘0, x) will bring about a smaller decrease in 
the norm. However, such alternate values of z may still be of interest. 

The choice of 


6= 2°'(8 — y) (27) 
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recommends itself for this purpose, for the following reasons: 

a. If A is hermitian, then y = — 8, and @ = 8. In particular, if A is symmetric, then 
8 = 0 or z, and, consequently, also 6 = 0 or x. Thus z = RF or — R, and T remains 
free of imaginaries. 

b. A is skew hermitian. Now y = r — 6, 6 = 8B — 2/2. For skew symmetric matrices, 
then, depending on whether 6 = 0 or x, 6 = — 2/2 or 2/2, so that z = — Ri or Ri, as 
desired. 

ce. Above choice of @ simplifies the calculation of the cubic somewhat. As may be 
inferred from Eqs. (23), (24), (25), the arguments 6 + y — 6,0+ 7 — a, 0 — B+ 4, 
@— 8+ a,26 —8+7, become, respectively, 2° (@ + y) —-6=¢9¢,2' (+ y7) - 
a=y, —¢, — vy, and 0. 

As the number of iterations increases, the super-diagonal elements in the pth trans- 
formed matrix A, tend toward zero. Consequently also b — 0, a, — 0, a,; — O for 
i <k <j. Simultaneously, F(R, 6) — RG(R, 6), where 


G(R, 0) = !c |’ R’ —2/e|[\|d| cos(@+y — 8) F (28) 
— a|cos(@+y—a@)R+|d—al?+ > (+4). 
3 


But the discriminant of G(R, 6) = O equals (| d | cos (@ + y — 6) — | a 
cos (06 ++ y — a))* —|d—al’ — ; , and is thus negative. Therefore, as the iterations 
proceed the largest real root R tends to zero, so that the matrix tends to its triangularized 
form. 

5. Computational procedure. The foregoing discussion indicates how the procedure 
may be applied in practice. 

A method of selecting the pivot 6 must be decided upon. Once the triangle of lower 
norm has been determined one may take as pivot the element of largest absolute value, 
any above-average element, proceed cyclically to the ‘‘next’’ element in that triangle, 
or take b in some other fashion. Pivots lying in the super-diagonal (j = i + 1) simplify 
the computation somewhat in that they cause the sums >>, in Eqs. (24) and (25) to 
vanish. With more experience better choices of pivots may become available. 

Next, the selection of angles @; of rotation must be made for insertion into the co- 
efficients of the cubic, and ultimately, for the transformation itself. Here, again, more 
experience will provide better clues. 

Thus, one starts with the calculation of N’ = p rr a;; |", N” = + er a;,; |", and 
the determination of the “upper” triangle by means of N’”’ = min (N’, N”’). For checking 


purposes, 


NA) = N’+N"+D, with D= > |a,; |’, 


as well as tr(A) = >°, a,; should be computed. 
1. Once the pivot b in A, has been selected, | a|,|b|,|e¢|,!d@/, | ai |, | aj |, | aes |, 


. . . . —1 . 
a,; | fori < k <j are calculated. From a = a’ + ta”, cosa = a’ | a} ,sina = a” 
a|', +++, and similarly, cos a,; , sin a; are obtained. 
2. Knowing @, ,7 = 1, 2, --- , m, there are computed the cos functions occuring in 


the C, . Thus, for example, 


cos (0, + y — 6) = (Je | | d |)~"[eos 0,(c’d’ + c’’d’’) — sin 0,(c’’d’ — c’d’’)). (29) 


t i 
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3. The coefficients C, of the cubic (21) are now determined, and the largest real 
toot R; of F(R; , @;) = 0 is found. This step is repeated for the other 6; chosen. Then 
the FR of largest absolute value among the R; is ascertained, as well as its associated 0. 

4. Now z, 7, tare computed by Eqs. (7), (17), and the elements of A,,, are determined 
by formulas (8)-(15). They are checked by Eqs. (16). 

la. A pivot is selected in A,,, , and the next iteration is begun, as described in step 1. 

As may be inferred from this computational outline, the procedure lends itself 
readily to easy programming on modern high speed computing machinery. 

A count of arithmetical operations reveals that approximately 8n + 8(7 — 7) + 8 
multiplications of complex numbers are involved per iteration, aside from the solution 
of the cubic. For an average value of (n + 1)/3 for 7 — 7 this amounts to about (32/3) 
(n + 1) multiplications per iteration. For each eligible element to act once as a pivot— 
which is the least one might expect—n(n — 1)/2 iterations would be required. This would 
result in » — 1 transformations of each diagonal element, at a cost of about (16/3)n* 
multiplications. This estimate seems reasonable in view of counts of (4/3)n* multipli- 
cations for the reduction of real symmetric matrices [2]. 

6. Example. The following simple example illustrates some of the salient features 
of the method. Let us take 





1 0 -2] 
A=/2 -l 2 (30) 
l2 1 oJ 
The characteristic values of A are 
M=—-2, wMH14+2, A= 1— 2. 
We note that tr(A) = 0, N?(A) = 19. Further, 
n+ Lief os, D= Y|\ [= 14. 

iat = 


In carrying out the reduction of A to triangular form the element of largest absolute 
value was taken in general as the pivot, and the angle of rotation @ was determined 


according to Eq. (27). 
Since N’ = 8, N” = 9, we have N’” = N’; the matrix A, as it appears in Eq. (30), 


has already its “upper” half in the proper place. 


1. The element b = — 2 is the pivot, whence a = 1, c = 2,d = 0. Further, a = y = 
6= 0,8 = 7. 

2. Next 6 = 27' (8 — y) = 2/2, whence gy = 2" (8 + y)'— 6 = 2/2, y = 2” 
(6+ 7) —a= 2/2. 

3. Therefore, C, = 4, C, = 0,C, = — 2, Cy) = 0, and F(R) = 4R* — 2R. Thus 


R = 277 6 = 2/2. | 
4. It follows that z = 27’ i, r = (2/3)'”, t = (2/3)'”7,.and 
[ 66667 + 1.885617 .57735 66667 + .47140i 
A, =| 1.15470 + 1.63299) —1 —1.15470 — 1.63299: |. (31) 


— .66667 — .471407 .816502 33303 — 1.885617_| 


| 
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la. Now the element ai? becomes the pivot. This leads to 6 = 72 22’, and F,(R) = 
.66667R® + 3.10395 R? + 2.06733 R — 7.36853. It turns out that R = 1.13916, and 


consequently, 


.66667 + 1.88561: 61638 + .71724i —.00428 — .32494¢ 
A, =| .75643 — .75390i —1.57773 — .64200i —.69069 — .50386i |. (32) 
L—.70462 — 1.74055 «1.53811 + .97889i 91106 — 1.24361: 
After eight iterations there is obtained 
95014 + 2.119921 — .22712 + .15730i .13567 + .69079% | 
A, = | -1.31007 — 280451 —1.97139 — .073421 24383 — 524671 |. (33) 
— .76294 — .14967i — .91102 + .58694¢ 1.02125 — 2.04650: 


The real parts of the characteristic values are thus already determined to better than 
5 per cent, while the imaginary parts are off by less than 6 per cent. 

7. Conclusion. A new method has been presented for the determination of the 
characteristic values of arbitrary matrices having complex elements. The example 
shown above, as well as additional computational evidence, justify the belief that the 
approach is a promising one. 


REFERENCES 


1. C. C. MacDuffee, The theory of matrices, Berlin, 1933 
2. W. Givens, A method of computing eigenvalues and eigenvectors . . . Nat. Bur. Standards, Appl. Math. 
Ser., No. 29, 1953 





276 BOOK REVIEWS [Vol. XIV, No. 3 


BOOK REVIEWS 
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Some basic problems of the mathematical theory of elasticity. By N. I. Muskhelishvili. 
Third, revised and augmented edition, Moscow and Leningrad, 1949. Translated 
from the Russian by J. R. M. Radok. P. Noordhoff, Ltd., Groningen—Holland, 
1953. xxxi + 704 pp. $10.05. 


The first Russian edition of this work appeared in 1933. Though the essentials of the author’s method 
of solving two-dimensional elastostatic problems have been repeatedly presented in the Western literature, 
no comprehensive treatment of this method has so far been available in any but the Russian language. 
While the number of Western scientists capable of working their way through a scientific note or short 
paper in Russian is fortunately increasing, a work of over 700 pages on a highly specialized subject is 
practically inaccessible to almost all. The Translator and the Publishers have therefore rendered a sig- 
nificant service to science by making this material available to Western scholars. 

The space available for this review does not permit detailed comments on the contents of the work 
or even reproduction of the table of contents which fills over ten pages. The following headings of the 
major parts (with summarized chapter headings added in parentheses) must therefore suffice to give an 
idea of the contents. I—Fundamental equations of the mechanics of an elastic body (analysis of stress 
and strain; basic equations of elasticity). II—General formulae of the plane theory of elasticity (basic 
equations and complex representation of solution; multi-valued displacements; thermal stresses). 
I1I—Solution of several problems of the plane theory of elasticity by means of power series (Fourier 
series; regions bounded by a circle; circular ring). 1VY—On Cauchy integrals. V—Application of Cauchy 
integrals to the solution of boundary problems in plane elasticity (regions bounded by a single contour, 
half-plane and other semi-infinite regions). VI—Solution of the boundary problem of the plane-theory 
of elasticity by reduction to the problem of linear relationship (sectionally holomorphic functions, half- 
plane and plane with straight and circular cuts). VII—Extension, torsion, and bending of homogeneous 
and compound bars (problem of Saint Venant, bars consisting of different materials). 


W. PRAGER 


Numerical mathematical analysis. By James B. Scarborough. The Johns Hopkins Press, 
Baltimore, and Oxford University Press, London, 1955. xix + 554 pp. $6.00. 


This book is the third edition of a well-known elementary introduction to numerical analysis first 
published in 1930. The treatment presupposes only calculus and relies heavily on examples worked in the 
text. The most glaring omission in the second (1950) edition has been remedied by the inclusion of a 
chapter on the numerical solution of simultaneous linear equations in which several of the best methods 
of solving such systems are treated and illustrated by examples. A further addition to this edition is an 
article on errors in determinants. Only the barest reference is made to the usefulness of methods for 
automatic digital computers. There are chapters on accuracy (I), interpolation (II-VI), differentiation 
and integration (VII), quadrature (VIII), algebraic and transcendental equations (IX-X), ordinary (XI) 
and partial (XII) differential equations, integral equations (XIII), normal error law (XIV), precision of 
measurements (XV), empirical formulas (XVI), harmonic analysis (XVII) and simultaneous linear 
equations (XVIII). 

The book is extremely well written and should be particularly useful because of its many examples 
and very valuable to students with only a moderate mathematical background. 

WALTER FREIBERGER 
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THE BENDING OF A WORK-HARDENING CIRCULAR PLATE 
BY A UNIFORM TRANSVERSE LOAD* 


BY 
WILLIAM E. BOYCE 


Brown University 


Abstract. This paper contains an analysis of the bending moments and deflection 
of a work-hardening circular plate under the action of a uniformly distributed transverse 
load. A segment-wise linear yield condition and the associated flow rules are used in 
order to avoid the unsound features of total stress-strain laws while retaining much of 
their mathematical simplicity. 

1. Introduction. For several years investigators in the field of work-hardening 
plastic solids have been aware of the physical shortcomings of total stress-strain laws 
(see, for example, [1]**). Consequently all recent theoretical work has involved in- 
cremental laws. However, since total laws are mathematically more tractable than 
incremental laws, solutions to many specific problems are still based on the total laws. 
In such papers a plausibility argument is given to indicate that the results are approxi- 
mately correct, despite their questionable derivation. 

W. Prager [2] has recently suggested a compromise which provides, to a large extent 
at least, the mathematical advantages of total stress-strain laws, while retaining the 
fundamental character of incremental laws. This is accomplished through the use of a 
segmentwise linear yield condition together with the associated flow rules. In order to 
make this paper self-contained a brief summary of this method is given in See. 2. 

In See. 3 the bending of a uniformly loaded, simply supported circular plate is studied. 
For uniaxial stress the mechanical behavior of the rigid work-hardening plate material 

0 


Pare 





A PR 7 


b— 





ie 


Fig. 1. 


is indicated in Fig. 1. For plane stress with the principal stresses ¢, and a, , the yield 
condition is assumed to be given by the polygon in Fig. 3. 

2. Review. Figure 2 depicts a simple model which represents the mechanical behavior 
of work-hardening plastic solids in plane stress with fixed principal axes. Let o, , 2 , 


*Received Oct. 24, 1955. 


**Numbers in square brackets refer to the Bibliography at the end of the paper. 
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€, » €, denote the principal stresses and strains respectively. The elliptical frame F is 
mounted so that it can undergo any translation in its plane. The perfectly smooth pin 
P is free to move inside the frame and may impart motion to the frame only when in 
contact with it. If there is no friction between frame and pin, such motion will always 
have the direction of the exterior normal to the frame at the point of contact. As the 
frame is displaced by the moving pin, the vector OP has the components co, and a; , 
while the vector 0C has the components ce, and ce, , c being the slope of the stress- 
strain diagram in Fig. 1. The elliptical frame corresponds to the yield condition of v. 
Mises, which in plane stress takes the form 


2 2 2 
01 + 62 — 0102 = %, (1) 


where a> is the yield stress in simple tension. 

A model of the sort just described clearly represents an incremental stress-strain 
law in which the state of strain at any time depends not only on the state of stress at 
that time but also on the loading history. It is precisely at this point that incremental 
laws may become extremely cumbersome, as the loading process must be followed 
step by step in order to determine the final states of stress and strain. 


The mechanical behavior indicated by the model of Fig. 2 will not change sub- 


CE, Oe 











Fia. 2. 


stantially if the elliptical frame is replaced by a polygonal frame, as in Fig. 3, provided 
that the latter has a sufficient number of sides. Either inscribed or circumscribed poly- 
gonal approximations may be employed. 

There is, however, a significant mathematical difference between the nature of the 
flows corresponding to the two types of frames. As long as the pin pushes against a given 
side of the polygonal frame, the translation of the frame is in the direction of the exterior 
normal. Consequently, the instantaneous position of the frame depends only on: (1) 
the time at which the pin made contact with the side in question, and (2) the instan- 
taneous position of the pin. Similarly, as long as the pin engages a given corner of the 
frame, the displacement of the frame coincides with that of the pin. Hence the instan- 
taneous position of the frame depends only on: (1) the time at which the pin first engaged 
the given corner, and (2) the instantaneous position of the pin. 

In this way a segmentwise linear yield condition and the associated flow rules result 
in the partial independence of the final strain from the loading history. This permits 
the use of total stress-strain laws “in the small’. Since the pin will usually contact 
































BENDING OF A WORK-HARDENING CIRCULAR PLATE 279 





1956] 


several sides and corners of the frame during the loading process, several types of total 
laws may be required in succession. Thus the fundamental nature of incremental laws 
is retained “in the large’. 

3. Bending of work-hardening plates. Consider a circular plate of radius R and 
uniform thickness h made of a rigid, work-hardening material. Let the plate be simply 
supported around the entire edge, and suppose that it is subjected to a uniform transverse 
load of intensity p per unit area. The usual approximations of the technical theory of 
the bending of plates allow the state of stress and strain at any point to be specified by 
the radial and circumferential bending moments, M and N, and the radial and circumfer- 
ential curvatures, x and \. The bending moments must satisfy the equation of equilibrium 


(rM)’ — N = —3pr’, (2) 


where the prime denotes differentiation with respect to r. In terms of the deflection w 
of the middle surface of the plate, the curvatures are given by 


k= —w”, A= —- w’. (3) 


On the basis of these definitions x and \ must satisfy a compatibility equation, namely, 





(rh)’ = x. (4) 
O g(N,cd) 

-- (M,) —+ 
D * 

2g T 

m,) A&M) 

= 
(M,cx) 

as-aMy- 





Fic. 3. 


The relation between bending moments and curvatures is indicated by the pin- 
frame model in Fig. 3, where the labels in parentheses must be used. It will be useful to 
speak of plastic states represented by the interior of the side CD, for example, as regime 
CD, and similarly for the other sides and vertices of the yield polygon. 

The problem now falls naturally into two parts: first, the determination of the 
critical load p, , the bending moment distribution, and the velocity field at the onset of 
plastic flow; second, the determination of the bending moments and deflection at some 
later time, that is, for p > po. 

At the beginning of plastic flow no work-hardening has taken place, and the analysis 
may be carried out by well-known methods of the theory of rigid, perfectly plastic 
solids [3, 4]. 

3.1. Incipient plastic flow. At the onset of plastic flow, M = N = M, at the center 
of the plate, so that the center is in regime B (see Fig. 3). In the neighborhood of the 
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center it is easy to show from the equilibrium equation that regime BC applies. At 
some-radius r = p, regime C is reached, and for r > p regime CD occurs. Since M(R) = 0 
the edge of the plate is in regime D. The yield condition and flow rules in these regimes 


are: 
3 (3? —-1)M,.<M<M,, N=2M,—M, d«k=da, (5) 


9 9 
r>p:0< MS 3i73 (3? — 1)M, , N = ai7a M, , dx = 0, (6) 


where d\ and dx are the increments of curvature corresponding to an increment dp of 
the load. 

The yield conditions (5) and (6) can be substituted into the equilibrium equation 
to give two differential equations for M, one valid in 0 < r < p, the otherin p < r < R. 
Their integration is subject to the boundary conditions M(0) = M, , M(R) = 0, so that 


lu, - 5 por’, r<p (7) 
M=4 
2 =) 1 (2 4 / 
7s M1 -- —) + = poR*(— — =), a 
giz M ( r g Pol | — pe rap (8) 


The transition point p and critical load p, are obtained from the fact that 





M() = 303 (3? — 1)M, , (9) 

because at r = p the plate must be in regime C. Combination of (9) with (7) and (8) 
furnishes the equations 

# a( sr 7 ae (10) 

Re 125i -1) oF gt 37 ote _—— (11) 


in the unknowns p/FR and M,/(p,R°). The solutions are found to be 


p ry Ae 

- = 0.4! 9 
p = 0-431, (12) 
Do = 6.66 a (13) 


It now remains to be shown that there exists a velocity field for the incipient plastic 
flow that is compatible with the stress field just obtained. This is necessary to insure 
that the preceding statical analysis provides the actual value of p, and not merely a 
lower bound. If the dot superscript denotes differentiation with respect to time (which 
here may be identified with differentiation with respect to the monotonically increasing 
load intensity p), such a velocity field must satisfy the following relations, the first 
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two of which spring from the flow rules (5), (6): 


1 
or See ey, | Fae 
14 
yp’ =0 > @, (14) 
w = 0, r=R, 
w', w’’ continuous at r = p. 
As is readily verified, the following velocity field satisfies these conditions: 
. r 
he<elek <x» 
. 2epR — p |’ = FP) ; 
wp = p. p (1s) 
. R-er 
Wo R — p/2 ’ ‘> 2. 


3.2. Work-hardening effects. For a work-hardening material, plastic flow will con- 
tinue only under increasing loads. For p > py , the region of transition between regimes 
BC and CD enlarges from a circle into an annulus which is in regime C. The situation 
is as indicated in Fig. 4. Zones 1 and 4 contain materia] which has been in regimes BC 


oO P R r 















- T(p) . ; ‘ 
w YY 4 } : 
! 1 ’ ’ : 
zone! oO Toy ® ‘ ® 
REGIME | ‘©; © ; 
( Fic. 4. 


and CD respectively since the beginning of plastic flow; on the other hand, the material 
in zones 2 and 3 is now in regime C but was originally in regimes BC and CD, respectively. 
Each regime must be considered separately. 

Zones 2 and 3. Elements of the plate in zones 2 and 3 are now in regime C but at the 
onset of plastic flow were in either regime BC or regime CD. Under many loading programs 
the relation between the total curvatures and bending moments for such an element 
is the same as if the element had always been in regime C, provided that the initial 
conditions at the appropriate point of BC or CD are replaced by those at C. Permissible 
types of loading are those for which the pin always remains engaged in a given corner 
of the frame, once it has become so engaged. After a solution has been obtained on this 
basis, it is necessary to verify that this restriction is indeed satisfied. Figure 5 will help 
to clarify this point. 

Since in regime C the pin is engaged in the corner C of the frame, the motions of pin 
and frame are identical. Hence 


dM = cdk, dN =cAdy. (16) 
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In view of the above observation, (16) may be integrated subject to the initial conditions 





ie - 
M=sn(38'"-1)M., N=3iMo, (17) 
giving 
| 2 
Me a7 (3/2 — 1)M, + cx, N = ai? M,+cx. (18) 
Nica 
M 
cK 








I. Initial position of frame. 
II. Position of frame when pin P engages corner B. 
III. Later‘position of frame. 


Eliminating x and \ among (18) and the compatibility equation (4) furnishes 


‘N’=M-N+35(2- 3M, (19) 
which must be solved simultaneously with the equilibrium equation (2). 

Boundary conditions are obtained from consideration of the transition points between 
regime C and regimes BC and CD. The boundary between regimes BC and C (zones 1 
and 2) will occur, say atr = o(p), when the pin reaches the straight line through C normal 
to BC; that is, when 


9 4 
M-N= — 317 (2 —3'”)M,, r=¢. (20) 
Similarly, the boundary between regimes C and CD (zones 3 and 4) is reached, at a 


point r = 7(p), when the pin touches the straight line through C normal to CD; that is, 
when 


2 


Now subtracting (19) from (2) and integrating with respect to (20) gives 


M 3(37-—1)M., r=r. (21) 


M-N= -; pr — 51 (3 — 3”)My | (22) 


1 wit.» ] 
73 (2 — 3'*)M, + a E po ~ 317 
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Substituting (22) into (2) and integrating subject to (21) furnishes 


2 " 1 is r o a” 
M = ai7a [3'? — 1]Mo + 503173) (2 — 3'”’)M,| 2 — “3 
. " . (23) 
+a [ara + 5-5] 
“ss *f FF 
and hence 
_ A as wal r " -§-§], 
N = M, + 2(3” 5 (2 3 ya 2 In 2 In +3 42. S| + 16 BP | 3: r’ r (24) 


Equations (23) and (24) can now be used to show that the condition on the type of 
loading mentioned earlier in this section is, in fact, satisfied. 

Zone 1. Elements in zone 1 have been in regime BC since the onset of plastic flow. 
The flow rule (5) can thus be integrated to give x = \, and then the compatibility equa- 
tion becomes r\’ = 0. Hence 


k= = f(p). (25) 
A further relation valid for regime BC is obtained from the fact that, as long as the 


pin pushes against the side BC of the frame, the displacement of the frame equals the 
component of the pin displacement that is normal to BC. This kinematical relation is 


expressed by 
cxAM + c\AN = c'(K’ + 2’) 
or, by writing out AM and AN in full and using (25), by 


M(p) — M(p.) + N(p) — N(po) = 2cf(p). (26) 
Now at the beginning of flow, (5) states that 
M(po.) + N(po) = 2M, ; 
therefore (26) becomes 
M + N — 2M, = 2cf(p). (27) 


The function f(p) is obtained by using the values of M and N at r = o as furnished 
by (23) and (24). Thus 


9 


] 
2cf(p) = 31 3(2— 4 )M,(2 In 2 In © = + $) - 7 (2 — 3 *)\M 


(28) 
? : (s ? — 267 — ) 
3 P T o 7) 
and 
y 2 
M+N= (3 _ soa)Mo + z73 (2 — 3” »M(2 = 4 S) 
3 5 * 2 . 
(29) 


1 4 
+ holar — 26 - § 
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Finally, substituting from (29) into the equilibrium equation, and integrating under 


the restriction that M be finite at r = 0 provides 


9 , 
M = 1(s _ 2,)M + mae (2 — 3' »M(2 in = + “) 
2 \ ee 2(: ee. (30) 


l , ; $) ] A | 
. ¢ P cee Se eee .~ 
+ 16 r( 3+ 20 3 g Pr. 
Then 


“(2 — 3' »u,(2 In = + A 
- =e (31) 


4 


l $5) ] ' 
9,2 _ 9,2 ; 2 
-+- 16 r( 3: 20 3 + g Pr 


It is clear that the symmetry condition WM = N at r = 0 is satisfied. 
Furthermore, (20) supplies a relation between a, 7, and p, provided that values of M 


and N from (30) and (31) are used; this equation is 


gc Pp Po 
, = 5 : (32 
R° R° P 


when (10) has also been employed. 
Zone 4. The material in zone 4 has been in regime CD since plastic flow began; hence 
ig x = 0. The compatibility equation then 


the flow rule (6) may be integrated, giving 


becomes 
ra’ +A = 0, 


so that 
I ” 
A = — g(p). (33) 
r 


The kinematical relation valid in zone 4 is obtained from the fact that the displace- 
ment of the frame equals the vertical component of the pin displacement, that is, 
dN =cad,y. 
The initial value of the circumferential bending moment is given by (6) while the 
circumferential curvature vanishes when the pin first makes contact with the frame. 


Thus in zone 4 


9 9 . 
N = sia Mo + ed = aia Mo + < g(p)- (34) 
3 3 r 
The function g(p) may be obtained from the value of N at r = 7 as given by (24). 
Hence 
Cc l “ e1/2 o l “ o1/2 ] 2 i = 
— gp) = ai72 (2—3°°)M, < fine 3172 (2-—3°°)M, + g P i 2} (35) 
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and (34) then gives 


2 ] 
172 M, pas 317? 


: oe ae . l p : 

N= (2 — 3°) =2 (8 — o) + - 2 (s' — o*. (36) 
3 TT 8 7r 

The equilibrium equation now provides an expression for M, in which the integration 

is carried out subject to M = 0 at r = R, as required by the conditions of simple support. 


Thus 


) ? . 
M = =a; NAC - =} - 37 (2 — 3”) ue (7? — o”) In 7 
\ / e T 


a ae ‘ 
+ go fr o ) In 


(37) 
f 


lp i 
R + 6r (k - 


A second relation between o, 7, and p is gotten from (21) by using the value of M(r) 


from (37), 
2M ' :) l ts o* ) T 
-(2 — 3'° — + - ; — —33] In 
3172 Dp? a 8 \ p? 32 
3° pk fi TR R (38) 


-at 2 — 37) Ma ( 5) r (8 . 2) 
giz (2 — 3°) 2 (1 — S) ng + G(T — ge) = © 


Now only the displacement of each element of the plate remains to be found. Using 
the definition of \ in terms of the displacement w, differential equations for w, one valid 
in each regime, are furnished by (25), (18), and (33) as follows: 





0 <r < oe: w’ = —rf(p) | 

. 9 
srs — w’ = N(p) — 3i72 Mor. (39) 
r <7 <= 2: w’ = —g(p) J 


Of course, the functions f(p), N(p), and g(p) are those given by (28), (24), and (35) 
respectively. The integration of (39) must be carried out subject to the boundary con- 
dition w = 0 at r = R and the requirement that w be continuous at r = o and r = +. 
The symmetry condition w’ = 0 at r = 0 is found to be automatically satisfied. The 


appropriate displacements are the following: 


cw l al/2 iz, =) o 
— (2 -—% >] In - 
MR 2(3'”") ' | Rk + Rk , T 


ro VR )~ of -§)| 

+(% e\24 9) aR \l- 2 (40) 
1 pk* p k o _ r = o a .@ 
AM pk oR fae th pt tee 


ae 1 (3 r - o ™ 1 Pa )| 
R? \2 R° R’ 2R°r° 
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o<r<r 
ow dg _ aia (2" =) r B(s a r) 
MR ~267)?— 3 | a id eae 7 
T 1 7° l or’ o” 
ee pe alee pee eee 2 
ie x Sox + 2 rR 7 =] (41) 
1 pR" p ! ee aes a 
+ 64 Mp po Lt Re tT RT SR BR 
a a a +4 
‘_ se ae a 
car <& &: 
i. > . -- 2 >2 > 4 4 
MR 3 T T/\R R 8 Mo po \r 7/\R R 


4. Numerical calculations and conclusions. Since the yield condition used in this 
investigation is represented by a polygon circumscribed to the Mises ellipse, the corre- 
sponding critical load p, at which plastic flow begins is an upper bound for the actual 
critical load predicted by the Mises theory. The latter quantity has been found by 
Hopkins and Wang [4] to be equal to 6.52 (M,/R*), whereas the analysis of Sec. 2 yields 
the present approximate value 6.66 (M,/R’). 

Equations (32) and (38) have been used to compute ¢/R and r/R for p/p, = 1.25, 
1.50. The results are shown in Table 1. Values of the normalized bending moments and 
deflection have also been computed for various values of r/R and are given in Table 2. 
These quantities are also portrayed by the solid curves in Figs. 6(a), (b), (c). The dotted 


TABLE 1 




















p/Po |i 1 1.25 1.5 
o/R 0.43103 0.38552 | 0.35194 
7/R 0.43103 0.63363 | 0.71801 





TABLE 2 




















r/R 0 0.2 0.4 0.6 0.8 1.0 

M/M, 1 0.96670 | 0.86679 | 0.68076 | 0.38850 | 0 

p/po = 1.00 N/M, 1 1.0333 1.1332 1.15470 | 1.15470 1.15470 
cw/M,R? 0 0 0 0 | 0 | 0 

— ' some ak" ~ | i 

M/M, 1.3179 1.2762 1.1511 0.90163 0.51128 | 0 

p/po = 1.25 N/Mo 1.3179 | 1.3595 | 1.4724 | 1.4310 | 1.3631 | 1.3214 
cw/M>R? 0.12223 | 0.11587 | 0.096802 0.066685) 0.033348) 0 
M/M, 1.6234 1.5735 1.4208 1.1089 0.62912 0 

p/po = 1.50 N/M, 1.6234 | 1.6734 | 1.7754 1.7132 | 1.5937 1.5059 
cw/M,R? | 0.24907 | 0.23661 | 0.19922 | 0.13984 | 0.070242) 0 


| | | | 
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curves are those obtained by Prager [2] with the Tresca hexagon inscribed to the Mises 
ellipse as the yield condition. The radial bending moment distribution is quite similar 
in the two cases, while the circumferential bending moment and displacement exhibit 
more substantial variations. 
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ON THE SOLUTION OF PROBLEMS OF 
DYNAMIC PLANE ELASTICITY* 


BY 
J. R. M. RADOK** 


Aeronautical Research Laboratories, Melbourne, Australia 


Application of complex variable methods to the dynamic equations of the plane 
theory of elasticity is shown to lead to a similarly lucid presentation of boundary value 
problems as in the case of the static equations. A stress function may be introduced, which 
is shown to be a solution of a fourth order partial differential equation which reduces 
to a generalized biharmonic equation of the type occurring in anisotropic elasticity, 
when it is required to find plane wave solutions. This stress function can then be ex- 
pressed in terms of two analytic functions of two complex variables and the stress and 
displacement components may be found in terms of these functions. 

Three problems are solved by this method to illustrate its advantages: the moving 
Griffith crack, punch and dislocations. 

1. Introduction. The introduction of the theory of functions of a complex variable 
into static plane elasticity of isotropic media [1] has considerably widened the scope for 
solving specific problems in this field. The reason for this development is to be found 
in the replacement of the indirect approach to the solution of boundary value problems 
by a direct one. In the present paper an attempt has been made to apply similar reason- 
ing to the dynamic equations of plane elasticity. 

The main difference between the methods of this report and those of [1] arises from 
the fact that the stress function here has to satisfy a generalized biharmonic equation 
which is of the same type as that occurring in the plane theory of elasticity of anisotropic 
media. 8S. G. Lekhnitzki [2] first introduced the theory of functions of a complex variable 
into this subject and showed that the stress function can be represented in terms of two 
analytic functions of one of two different complex variables each. 

Some of the results obtained below were deduced by different means in 1937 by 
E. C. Titchmarsh [3] when studying the theory of Hilbert transforms. 

2. The basic equations. In the case of plane elasticity the equations of motion are 


Oc, Ot, ru 
Ox Oy ot (2.1) 
OT x, e .. av 
Ox dy P at 


The stress-strain relations are the same as in the case of static elasticity 


, Ou Ov 
gg = hO + 2p — , o, = r\0+ 2-—, 
OX Oy 99 
ou Ov Ou Ov 
Tx, = A({ a ), Gm = 4 => 
Oy Ox Ox Oy 
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The first two of these may be rewritten in a form which will be found useful later on 


Ou r ov A 
a = — 3 2a— = — 
¥ Ox Cs 2(r + a 6% + % ‘ dy " 2(r - +p y (% +o), 
(2.3) 
§= (c, a 
20 "9 yy wet % 
Finally, the equation of compatibility must be satisfied 
3 ri \, 2(A + p) ri (2 av) 
5+ <5) »- ar = (S+ =) =0; 2.4 
‘e F 5ypfitt + — \ pon oF lar t ay * G4 
using (2.3), this may be rewritten 
[ 9? a 3 
: ) See Teme f = (Q A? 
lax? + ay” r + Qu at 1oz + Cy} 0. (2.4 ) 


Differentiating (2.1) with respect to x and y respectively and adding one finds by (2.2 
and (2.3) 


oe “Oz d°o, 3 ’ n 
eli i Pies e ase (oc, me Ty); (2.5) 
Ox” OY Zu Ot 
ae f x o” po 
a a, Os oy) (: 3 —— oS Lan = 0. , 
Ox OY t oy) + Ou” T oy =pot] * (2.6) 


The complex representation of a certain type of solution of Eqs. (2.4’)—(2.6) is the 
object of the next section. 

3. The stress function and its complex representation. Equation (2.5) may be 
rewritten 


{ 9? _ Pp é\ is a 2 a*| . (3.1) 


Ane 4.2 ajz(o 
Ox” 2u ot 7" lay" 2u ot) 


and it is seen that this equation will be satisfied by a stress function U defined in the 


42 42 FY a2 ; 
0; * {a me oe \, : Co, = f é a & : hw. (3.2) 


following manner: 


lay" Qu at’) lax Qu ot 
Using (3.2) one has 
a 3° p id 
me fas ass — = ofl (3.3 
ten laat ag fae ) 
and (2.6) and (2.4’) become 
{a ——*- a ty +=} = 0, (3.4) 
(Ox oy pol ) Ox Oy) 
0° a p a’) (a? 4 ™ 
——— + = ta pee _ = 0. (3.5 
‘2 oy” A + 2y Ot") \Ox t ay im 0 hy ) 


It is thus seen that U must be the solution of a generalized biharmonic equation, 
formed by the product of two differential wave operators. The sum of the shear stress 
7., and the mixed second derivative 0°U/dx dy must also satisfy a wave equation. It is 
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clear that in the case when the time ¢ is not a variable, the function U becomes Airy’s 


stress function. 
In what follows attention will be restricted to problems for which the following 


transformation holds true 

§é=2z2+ct, n= Y;3 (3.6) 
such a transformation corresponds to the occurrence of disturbances moving with 
velocity ¢ parallel to the axis Ox. It is obvious that a more general approach could be 
adopted by use of the transformation 


é=r7+¢,l, n=yxCl. 


However, this has not been done here because it tends to complicate the following 
equations without introducing any basically new aspects; in particular, a rotation of the 
system of axes about the origin would lead to the same results. Using (3.6) Eq. (3.5) 


becomes 





# (4 Mas) 4 Ay ' 
‘& (1 A+ =) + oe . m ” an vila (3.7) 
with the real solution 
U = F,(z,) + File.) + FiG,) + F.@,), (3.8) 
where 
2=é&+ $17, Z=€& + 827, (3.9) 
s, being the roots of the characteristic equation of (3.7) 
2 2 
pe 2 pc 2 
a ins Es - 1 
{1 + eh 4 #| 0, (3.10) 
i. 
2 1/2 2\1/2 
il pc aa _ i pc" oe 
ss; = i{1 —_ d we =) —_ 1B; ’ So ed (1 — =") —_ 1B. . (3.11) 


The functions F;(z;) are analytic functions of the complex variables z,; and F,(Z,) are 
the corresponding conjugate functions. 

It is now possible to express stresses and displacements in terms of derivatives of the 
functions F;(z;) and hence to reduce boundary value problems of dynamic elasticity to 
problems of complex function theory. 

4. Stresses, displacements and boundary value problems. In order to simplify the 
later results, introduce the notation 





dF, dg dF, 
> an ele on ie - ; (4.1 
. dz,’ ? dz,’ v dz,’ ws ) 
Then, by (3.6), (3.8) and (3.9), 

aU 
oy = 2Rele’ + v'}, 
PU ° 9 2 
on = 2 Relsig’ + sy’) = —2 RelBiy’ + Biv’), (4.2) 
a°U 2 0°U 


ae OF’ 
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while (3.2) take the form 


. 4? 2 yr 
-(1-2)Z0-la+afu, 
Zu 0& a dé 7 
(4.3) 


f 0° p C or \ 


\dn” 2u OF, 


o, = = 


Hence one finds the following expressions for the stresses and displacements: 


o, = (1+ 8) Rey’ + VY), 


o, = —2 Re{[8i + (1 — B)]e’ + [1 + Bl¥’}, 
go, +o, = —2Re{(B; — B)e’}, 
J. 74 10+8% , | , 
T,, = —2 Re 1) Bip + = 8, ae fo (4.4) 
, we 
= a _2 Pate oy + B2)¥/, 
Bh | 2 ) 
mwJ..,10+8,| 
y= = 0) Bue 5" 8, vfs 


where the last two have been obtained from (2.3) using the first three and integrations, 
while the expression for the shear stress r,, follows from the displacements by means 
of differentiations. In this work certain constants and linear terms have been omitted, 
since they only refer to rigid body motion. 

It is readily seen that the formulae (4.4) satisfy the differential equations of Sec. 2. 
They have been deduced previously by Titchmarsh [3] when studying Hilbert trans- 
forms, using entirely different methods. Nevertheless, they have not, to the author’s 
knowledge, been used for the reduction of boundary value problems of dynamic elasticity 
to problems of complex function theory. 

Certain types of boundary value problems of dynamic elasticity have been discussed 
in detail in Chap. 3 of Ref. [1], and it is shown there that the solutions are unique. In 
general, boundary value problems will be stated in physical terms and their reduction 
to function theory by means of the formulae will present no difficulties. As in the static 
theory of isotropic media full use can be made of conformal mapping, except that, just 
as in anisotropic theory, the occurrence of the two variables z; requires the use of two 
mapping functions. In fact, it is obvious that a boundary L in the z-plane, the original 
plane, will be distorted in the z, and z, planes. The statement of the boundary value 
problems in terms of two functions of z, and 2, respectively therefore will require conformal 
transformations from the z, and z, planes on to the circle y in the ¢ plane, where only 
one variable occurs. 

These remarks will be further elucidated in the next section which deals with a 
fundamental problem of the type J (stresses given on the boundary)—the problem of 
a finite length crack moving through a solid. This problem, solved by different means 
by E. Yoffe [4], is of some interest in micro-elastic work. Another problem, readily 
treated by this method, is that of a parabolic punch moving along the boundary of a half- 
plane. This problem is solved in Sec. 6. The problem of moving dislocations is briefly 


discussed in Sec. 7. i 
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5. The moving Griffith crack. Consider the case where uniform tension p acts at 
infinity at right angles to the direction of the slot, extending in the £, 7 plane from 
— ato + a. The boundary conditions for this problem have the form 


= 0 
~ E| <a 
f, = 0 (5.1) 
c, =P, cs, = 0. 
Using (4.4), one finds from these conditions 
Rele'(z:) + ’(@.)} = 0, 
(5.2) 


-2 Re fate) + L+H yoeyh = 0, 


and the problem is to find analytic functions ¢(z,), ¥(z.) which give the required stresses 


at infinity and satisfy these boundary conditions. 
Obviously the uniform loading at infinity will be caused by linear terms in ¢ and y; 


therefore one may write 
g(z:) = (B* + C*)z, + gla), ¥(z2) = (B*’ + C*’)z2 + Poles) (5.3) 
and express the constants B*, C*, B*’, C*’ in terms of the loading at infinity. By (4.4), 


one finds 








2) 2) (@) (a) (a) 
OR is aes +a, OP: ge ee o, + Ss... 5.4 
E oe —s)’ » ~T+Ht oe —s ’ aa 
C* = = , C=. 


Substituting (5.3) in (5.2) the modified boundary conditions become, since C*’ = 0, 
Re{ei(z:) + Wi(z)} = —(B* + Bt’), 


f ef 
Re isBygiles) + AEP yiie)) = B,C*. 


(5.5) 


Since in the present problem the boundary is part of the real axis, it is the same in 


the z, 2, , 2. planes and the mapping function on to the unit circle is 


ee 1 | . 
—— Si _ e 6 
a 2 (: + 1) 6s 
so that 
at=2.—-(e.-—a°)'” (5.7) 


Normalizing and transforming (5.5) and using the well-known result of function theory 
that an analytic function whose boundary value is constant is a constant throughout 


its range of definition, one finds 


eS) + pos) = —a(B* + B*’), 


2 ¢ on ¥ (5.8) 
vs) = aBC*. 


(1+ 8: 


B2 


nn 





Ce ee 


aes + 
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Consider now the special case 





1. =P, of =ty =0, 
i.e., 
B= -s9 ms 8 -itatae (5.9) 
C*=0 
Then, from (5.8), 
ap(1 + 83) 48 ,B.pa 








eS) = —T py — 4p,’ O -Gyeya +B) — 484] (5.10) 


and, after integrating with respect to ¢ and using (5.7), 





. pi + 83) f 2 2\ 1/2) 
Go\t1) = —7 2\2 -_— (zi — a) Sy 
(1 + 63) 48,8. (5.11) 





fe. — (23 — a’)'”*}. 


1,(Zo) = ’{ rr 
oles) (1 + 62)[(1 + 82)” — 46,82] 
The expressions for ¢(z,) and ¥(z.) may now be written down using (5.3). The stresses 
and displacements then follow from (4.4). 

In order to check the solution, the stresses along the real axis will be calculated. By 
(4.4), the stress distribution is given by 





Pp a 2 2 “1 
= > 9 — | 4) B; — ——-— 5-5 
ie 1 + 62)" — 46,8, ” \! a ee al (2; — a’)’ | 
22 1 
sai 18,8 ] = (2? a 2 ‘ 
4 — a) (5.12) 
i 1 + B3)p mH 2 13,8 zo | 
si 1 + 6) — 46,8 Ke\(1 F Be) | To - 1+ he-a)y 
—2p(1+ 636, , J 24 2 | 
fan © aT ae a ee hs OS OO 
(1 + 62)" — 48,8. \(z; — a) (2. —a)’’) 
and hence along the real axis 
28? + 1 — 62)\(1 + 62) — 46,6 P for 1€| <4, 
= - ————— J, 
° (1 + 6) — 48,8, gE 
0 for || <a, (5.13) 
o, = 4 ‘ 
En 
gam for || >a, 
T, = 0 for allé. 


It is thus seen that the solution satisfies all the conditions of the problem. 
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Next consider the transition to the static problem. By (3.11), 


a ee ‘ 
a, = (1 | = Fr +2 + MC); 


2\ 1/2 2 
a, = (1- =) =1-5- +06), 


and therefore 





i (26; + 1 — 63)(1 + 83) — 48,8. _ 
1m = —l 





pen (1 + 62)" — 48,8; 
and (5.13) gives 
—p for |x| <a, 

= ‘i \ 

—of1 _ Go ay) for |x|>a, 

f 0 for |x| <a (5.14) 

v, = la aE in for |x|>a, 
Ty = 0 for all z, 


which is the solution found by C. E. Inglis in 1913. 
The expressions (5.13) agree fully with the results obtained by E. Yoffe [4] whose 
analysis and formulae are easily related to those of the present report when it is realized 


that in her notation 

Soni ig(z,)B i(1 + B3) 

fe) =< g@2) = ——~5—— ve), y=86., B = p,. 

rm 2u 

The physical interpretation of the functions ¢g, y is deduced from [4], which constructs 
them from simple wave solutions by means of Fourier integrals. It is thus seen that 
y, w are associated with longitudinal and shear waves respectively travelling through 
the medium with the velocities 


Xr 4+ Qu 1/2 mn 1/2 
Cy, = (A+ ’ Cc. = es . 
p p 


6. The parabolic punch moving over a half-plane. It will be assumed that no friction 
occurs. Let the profile of the rigid punch be given by 


2 
x - 
f(z) =d-—- OR’ R large, (6.1) 


and assume that it is in contact with the elastic half-plane over the segment | | < 1 of 
the axis. Let the remaining part of the real axis be unloaded. Under these conditions 
the quasi-static boundary value problem takes the form 


Tz = 0 for all points of the ¢ axis, 
0 forj/é|>J, (6.2) 
f(§) for|é| <1. 


Oy 


v 
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By (4.4) these conditions lead to 


7 1 ( B3)° re, 
3, {e'(z:) — ¢ (2,)} 4 a — {W'(z.) — W(z)} = 0 for allrealé, 
v'(z:) +e) + We) + WG) =0 for |¢|>1, (6.3) 
, 7 Il ] B; ? rs ° . - 
Biig'(a) —e¢@}+5 a 2 tw’'(z.) — pW (%)} = Quf’ for |é&| < J, 


where the last equation has been obtained from (4.4) by differentiation with respect 
to £. These boundary conditions must be understood in the sense that they hold true for 


from inside the upper half-plane n > 0. 
A solution of (6.3) will now be sought in the form 


g(z:) = Ax(z,), ¥(2.) = Bx(z2), (6.4) 
where the initially arbitrary constants A and B are to be determined in such a way that 


the last two conditions (6.3) lead to a simple problem of linear relationship [1]. 
The first equation (6.3) gives then 


6 i 2\2 ; —/; 
Ee + ag, “! + 2) |e —x]=0 
for all real £, and hence 


B = 
AB, + 13 (1 + 63) = 0. (6.5) 


Similarily the last two conditions of (6.3) lead to 


(A + BY(x’ +x’) =0 for {é/ >, (6.6) 

and 
8, +1 + 6 Ik’ — x] = Quif’ for |e] <1 (6.7) 
AB, + 55 2) Ix’ — x’] = 2ut, < I. 4 


The constants A, B may now be determined from (6.5) and the condition 


3 , ” 
AB, + = (1 + Bi) =1 (6.8) 


deduced from (6.7). One thus finds 


1+ p; 13, ' 
A=-—",, B= (6.9) 
B,C — a) | — Be 
Note that in the presence of friction (6.5) would no longer be valid. 
It follows from (6.6) that the function x will be holomorphic in the entire plane cut 
along | | < J, if it is extended into the lower half-plane by 


— 


x(¢) = —x(¥). (6.10) 
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Then, as ¢ — ¢ from 7 > 0, f~ § from 7 < 0 and 
XH>xO™, xO>-x&. (6.11) 


Hence, using (6.8), (6.1) and (6.11), the boundary conditions (6.6) and (6.7) take the 
form 


it (6.12) 
x +x = —S for || <i. 


Using the Plemelj formulae [1] and taking into account that in view of the smoothness 
of the punch profile no singularities will occur, one has 


x(t) = —u(P — 2)? [ t dt 
Sr pa 2 2) 1/2 
Rr - (¢— $l ) (6.13) 








=F bg py”? — g}. 


The solution of the problem follows now from (6.4) and (6.9): 








, \ (1 1 2 2\1/% 
ve) = EB RB Ia a1 
vlc (6.14) 
ln Ye op _ de 
v'(z.) = — BS : (2 I’) ze}. 
The stresses are given by (4.4) and will here only be written down for the £ axis: 
= - Aa Be - An” for |&| <1, 
, (6.15) 
{48,8. — (1 + B3)(1 + 26; — B)} mu gy¥2 
7 ee f < 1. 
OC; (i — 2) LE ‘) or |—E| <1 
The total force acting on the punch is 
oe _ _ 148,8. — (1 + 63)*} _u Pa 
y=| oa = _ ma © (6.16) 


When the velocity c with which the punch moves in the direction & tends to zero, 
the above results reduce to those of the static theory, as was to be expected. 

7. Moving dislocation. In [5] the author looks for a solution of the differential 
equations of Sec. 2 which is to satisfy the following conditions: 

= 0 ; 
%y : for 7» =0 
[uJ], = const. 

where [ ], denotes the change undergone by the quantity inside the brackets for one 
circuit of L. It is easily seen from (4.4) that it follows from the first of these conditions 
that 


g(t) = —vé), 
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while the second is satisfied by the functions 


¢g = blogz, , y = —b logz,, 
which satisfy also the first condition. Hence, by the last two expressions of (4.4), 
Re 1 
u= am Hog = 9 (1 ote B>) log ad fy 
m 2 


| (7.1) 
Re . 1 (1 + 3)? } 
vy = —— 1b§B, logz, — = — — log ze 
cs © 2 Bs ass 9 
and, since 
log z = log (x + iy) = log (2* + y*)'”* + cartan y/z, 


the displacements (7.1) agree with those given by formulae (12) and (13) of [5]. 

In the work of this section, as also elsewhere in this paper, integration constants have 
been omitted. 

8. Conclusions. The introduction of the complex variable approach into the study 
of the dynamic equations of plane elasticity leads to a simplification of the statement 
as well as of the process of solution of problems dealing with surface waves. In fact, 
the relevant solution of the differential equations is found to be expressible in terms of 
two analytic functions, corresponding to the two types of waves, occurring in elastic 
media. As in static elasticity, the stress and displacement components are easily ex- 
pressed in terms of these two functions and in this way boundary value problems can 
be reduced to problems of complex function theory. 

The method has been applied to three problems some of which were previously 
solved by other more complicated methods: the moving Griffith crack [4], the moving 
punch and moving dislocations [5]. 

Not added in proof. Since completion of this paper the author has discovered that the 
problem of the moving punch has been solved by a similar method in [6]. 
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STEADY FLOW OF NON-NEWTONIAN FLUIDS 
THROUGH TUBES* 


BY 
A. E. GREEN anp R. 8. RIVLIN 


Brown University 


1. Introduction. The stress-rate of deformation relation for an isotropic fluid, in 
which the stress at any point is assumed to depend only on the velocity gradients at 
that point, has been considered by Reiner [1]** and by Rivlin [2]. If the fluid is in- 
compressible and the stress components are expressible as polynomials in the velocity 
gradients only, then the stress tensor can be expressed, using Cartesian tensors, in the 
form 


t;; = Od;; + Vdi,d,; + pd;; , (1.1) 


where d,; is the rate of deformation tensor defined in terms of the velocity vector v, by 
dj; = 40;,; + 9;.,); d;; = 0. (1.2) 


p is an arbitrary hydrostatic pressure and a comma denotes differentiation with respect 
to one of the coordinates x; . Also, 6 and WV are polynomials in the invariants IJ and III 
defined by 


II =3d,d;;, Il = 4didudy; . (1.3) 


tivlin [3] has discussed the steady rectilinear flow of a fluid, for which the stress is 
given by (1.1), through a uniform tube of circular cross-section and has shown that 
such a rectilinear flow can be maintained without the application of body forces, without 
making further assumptions regarding the form of the dependence of 6 and W on J] and 
III. Ericksent [4] has shown that this is not generally the case when the cross-section 
of the tube is non-circular and that, in such cases, the maintenance of steady rectilinear 
flow through the tube will require that, in addition to a uniform “pressure gradient”’ 
along the tube, an appropriate distribution of body forces be applied to the fluid. An 
exception occurs when WV = k®0, where k is a constant. In this case steady rectilinear 
flow of the fluid through the tube can be maintained by a uniform “pressure gradient” 
alone. 

In the present paper, a suitable distribution of body forces, for the maintenance of 
steady rectilinear flow through a non-circular tube, is calculated. This distribution has 
non-zero components only in the planes normal to the direction of flow and these com- 
ponents are independent of position along the tube. This result suggests that if a uniform 


*Received November 16, 1955. The results presented in this paper were obtained in the course of 
research sponsored by the Office of Ordnance Research, Department of the Army, under Contract No. 
DA-19-020-3487. 

**Numbers in square brackets refer to the Bibliography at the end of the paper. 

+We are grateful to Dr. J. L. Ericksen for discussing this result with us and for allowing us to see the 
typescript of his paper prior to publication. 
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“pressure gradient”’ is applied to the fluid, and no body forces are applied, a steady flow 
will be produced which consists of a rectilinear flow along the tube on which is superposed 
some flow distribution in planes perpendicular to the length of the tube, independent of 
position along the tube. This flow distribution is calculated for a tube having an elliptical 
1 fluid whose stress-rate of deformation relation has the form (1.1) 


cross-section, for ¢ 
a(c + II), where c is a constant and where a is a sufficiently 


with 6 constant and ¥ = 
small constant so that the governing equations determining the flow can be linearized 


with respect to a. This fluid is a nearly Newtonian fluid which departs from the ideal 
Newtonian law (obtained by taking a = 0) in a particular way. The choice of this 
particular form of departure is governed by considerations of mathematical tractability. 
It is found that for such a fluid, a uniform pressure gradient will produce in the elliptical 
tube a rectilinear flow on which is superposed four similar vortex-like flows in planes 
normal to the length of the tube—one in each quadrant of the elliptical cross-section. 
These have the same signs in diagonally opposite quadrants and opposite signs in 
adjacent quadrants. Although the quantitative results have been obtained for a par- 
ticular type of non-Newtonian fluid, we may expect them to be qualitatively valid 
much more generally. 

2. Rectilinear flow. Assuming zero body forces Ericksen [4] has examined steady 
flow for which the velocity components are 


vy, =v, = 0, v, = 2f(z, y), (2.1) 


in a fluid for which the stress-rate of deformation relation is given by (1.1), where co- 
ordinates x; are replaced by 2, y, z. Here we reconsider flow given by (2.1) but retain, at 
present, body forces. From (1.1) and (2.1) we have 


0, 0, f, 
d;; 0, 0, %. 3. 2.2 
> To 0 
is oie 0 
d;de; = || fefy; ae 0 |, (2.3) 
0, 0, fet | 


where suffixes denote partial differentiation. Also, from (1.1), (1.3), (2.2), (2.3) it follows 
that 
W=f.+f,, II =0, (2.4) 


ti, = Wf: +p, tie = Vf.f, , to = Vf, +p, 
(2.5) 


bh: = Of. ’ bs = Of, ’ t33 = Vf: + fo + P- 


If we denote body force components by (X, Y, Z) and the (constant) density by p, the 


usual equations of motion reduce to 
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fr) 


~~ ee oes . 
pX + ax + ax (Wf2) + ay (Wf.f,) = 0, 


pY + P+ SWS) + 5 Of (2.6) 


I 
= 


0 oP 


fe] 
Z+ P+ sf.) +5, (Of) =, 
if we use (2.5) and observe that inertia terms for the flow (2.1) vanish. 
We now seek for a solution in which Z is zero and X, Y depend only on z and y. It 
follows from the third of equations (2.6) that 


Pp wid Pox, y) + rZ, (2.7) 
where X is a constant such that 
5 (Of) + 5 Of.) +2 = (2.8) 


The function Y depends only on the invariant JJ which is given by (2.4). Hence, writing 


q=p +34 | vdll, (2.9) 
we have 
9g _ 9Po 
Ox = Ox + Wf fez + Ste, 
9g _ Po 
oy ors oy ' WS fer + Safer)- 


It then follows that the first two of Eqs. (2.6) may be rewritten in the form 


4 99 | oO a _ 
pX + ox + [2 (Wf,) + ay (Wf,) i — 0, 
(2.10) 
, , Og - 
pv + 24 | 2 Wf.) + a 2 ws |p, = 0. 


Ericksen [4] has shown that when X and Y are zero, Eqs. (2.10) imply that q is a 
function of f only, and, in general, they are incompatible with Eq. (2.8) for f. When 
X, Y are non-zero a possible solution of Eqs. (2.10) is 


pX -|2 wf) +2 5 BS) +: ‘aly. 


(2.11) 


Il 


pY 


-[200+ 20044). 


where g is now a function of f only. The flow given by (2.1), where f satisfies the differen- 
tial equation (2.8), together with a suitable boundary condition, can then be maintained 
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provided we apply body forces whose components (X, Y) along the (x, y) axes re- 
spectively, are given by (2.11). 

From (2.8) and (2.11) it is seen that the body forces vanish when 
(a) VW = kO, q = kXf, where k is a constant, or 
(b) fis a function of a linear combination of x and y, or 
(c) fis a function of (x? + y’)}, 
in agreement with Ericksen [4]. 

If the boundary of the cross-section of the tube is symmetrical about both the z- 
and y-axes, we see that equation (2.8) for f, together with the condition that f = 0 on 
the boundary, implies that f is an even function of x and an even function of y (6 depends 
only on f: + f,). Since dq/df in (2.11) is arbitrary we see that, in general, the flow may 
be maintained by a body force X along the z-axis which is an odd function of x and an 
even function of y; and by a body force Y which is an even function of z and an odd 
function of y. This suggests that rectilinear flow down such a tube may be possible 
without body forces, if in addition there is a flow in a perpendicular direction in each 
quadrant of the xy plane, the flow in any one quadrant being the mirror image of that 
in the immediately adjacent quadrants. In Sec. 5 this is confirmed for a special type 
of fluid flowing in a tube of elliptical cross-section. 

Suppose now that the fluid is such that the restriction (a) is satisfied approximately. 
We put 


Y=kO+ov’, q=hnrf, 2.12) 


and Eqs. (2.11) and (2.8) then yield 


i fa 0 1 
pX = —) 57 (w’f.) + ay | ‘ ”| Jas 
| . (2.13) 
oY = -a% (W’ f.) + = (W IDS ’ 


where W’ is a function of the invariant JJ, and a is a sufficiently small parameter. The 
value of g can be chosen arbitrarily and we choose to take the particular value in (2.12). 
In general, therefore, body forces (X, Y) of order @ are required to maintain the flow 
(2.1). For example, if the fluid is nearly Newtonian in such a way that 


0 = 2u, (2.14) 


then X, Y are O(a). 

In the next section we examine the modifications which are required in order to 
obtain steady flow in pipes, when body forces are zero, for a general fluid for which 
Eq. (1.1) holds. 

3. General steady flow in tubes. We seek a solution of the fundamental equations, 
with zero body forces, in which all the velocity components are independent of time ¢ 
and the coordinate z. We therefore introduce into (2.1) non-zero values of v, , v. . In view 
of the incompressibility condition d;; = 0 we may write 


Vv; = —y, 5] (3.1) 
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y and f being functions of z, y only. From (1.2) and (3.1) we have 
| ter Wat, fe 
dis = || Mee — Yue) cs & 
Se jf... OF 
DimAms = Wry + (Wee — Vo) +2» 
dindur = f.f,; 
DimOms = —feWey + 3h (Vez — Ww); 
don dms = Wry + UYez — Ww) + Si; 
domOms = SyWey + 2h (Vez — Vw)s 


3 mA m3 = fi + pb ’ 


and the strain invariants (1.3) are then 


TT = fo t+ fi t+ Viv t+ Yes — Vw)’ 
IIT = (fi — f2)Wev + $f bex — Vw) 

From (1.1), (3.2) and (3.3) we see that the components of stress are 
ty = —Ov., + Vf Vi, + vee — Ww)’ + fz} +2, 
tio = 30(vse2 — Ww) + VSS, ; 
too = Oey + Uden + Hee — Vw) + fi} +B, 
tis = Of. + Vl —Seber + 2h (Vez — Vw}; 
tos = Of, + Ul fuWer + 34S2(Ver — Vw)}; 


ts, = Vf +f) +p?. 
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(3.2) 


(3.3) 


(3.4) 


(3.5) 


Using (3.1) and (3.5) we may now write down the equations of motion. A more con- 
venient form of these equations may, however, be obtained by introducing the complex 


variables ¢, ¢* where 
f=z2+y, f=r7-y. 


Then, observing that 


(3.6) 


(3.7) 
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we find, from (3.5), that 


by ot. loo => 2p a AV (f fy + QU errWrere), 
tr — toe + ite = MOVs. + AV fp. , 
ths + Stes me 20f;. + 4iW f pProre ; 


3 =~ p+ 4Vf;fy. 


Also, using (3.1) and (3.6), the equations of motion become 


re] re] 

art (1, + too) + ag (ty ‘ins too T 2tti2) = 40( Vr Wrers a VreWere), 
Olas o P @ e = ° e = 

ny > of tis + tes) + ace (Lis — tts) = 4tp(frbs> — fro), 


so that with (3.8) we then have 


Op , , a — P , 
ort a 2 ar* WS fies aa QWerWrers) } 


oi. - mae 
2 {40 Y pepe + Uff} = 2o(Vedeere — Vrodees), 


o¢ ° 
Po? top. + Dif Vpn} 
= _ re ~ r rere 
Oz ' ” or t 73 d J re. ? 
D , beg oe 
+ 2 ro [Ofs — 2Wfredrr} = 4ip(frpy. — Srey). 
Of 


From these equations we see that 
| 
p=rze+ Pol f, 7). 


where X is a constant and 7, is a function of ¢, ¢*. Hence 


fe} ; o— a! ae : a oa 
+2 rr {Of rs + WWWfpWrere} +2 ae* (Ofs — 2WfreYee} 
= dip(frbre — fred), 
OD ee 
are 2 ae Wf irchre + WWeeWere) } 
>) a (2 ee, a ae 7 
+ 2 ae (VOW sere + VS re} = 2o(veveer — VreWrre). 


We also note that the invariants in (3.4) may now be written in the form 


IT = A(frh es + WerWeers), 


ITI = 4i(fjeYeers — f2eWe:). 
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(3.8) 


(3.9) 


(3.10) 


(3.11) 


(3.12) 


(3.13) 


(3.14) 


(3.15) 


Equations (3.13) and (3.14) are three equations for the three unknown functions 
] ] 


Po , f, ¥. In addition, if we impose the condition that the velocity is zero on some closed 
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curve G(¢, ¢*) = 0 in the (2, y) plane then 
f =0, v7; = 0 [G(¢, ¢*) = 0}. (3.16) 
In the next section we consider the solution of Eqs. (3.13) and (3.14), subject to boundary 
conditions (3.16), for a fluid for which 
VY = k0 + aV’. (3.17) 


This fluid was discussed in Sec. 2 where it was seen that, when a = 0, a solution of the 
problem was possible with y = 0. 

4. Solution by successive approximations. We suppose that W is given by (3.17). 
We seek a solution of Eqs. (3.13) and (3.14) in power series in a, based on the solution 


L = 0 when a = 0, assuming that a is a parameter which is sufficiently small, but we 


y 
only carry out detailed calculations for terms of 0(a). Thus, we assume that 


rapt, Jeartapt >, (4.1) 


so that, from (3.15) 


IT = 11, + all, +-:-, 





(4.2) 
III = alll, +---, 
where 
II, = Ah;h;. , IT, = A(hygye + hygs), 
(4.3) 
TTI, = 4i(hrgrere = hegre). 
Assuming that 0 and W’ are analytic functions of a, we then have 
OUT, 111) = @, +a0,+-:-, (4.4) 
WUT, 112) =V¥e+-::, 
where 
Q@ = O(II, , 0), Vv, = Vv'(II, , 0), 
(4.5) 
, _ zy 90(1T, , 0) aO(IIo , 0) 
O = 1, oll + HT; oll 


Hence, from (3.17) and (4.4), 
VY = kO, + alkO, + Wi) + --. (4.6) 


We now substitute (4.1), (4.4) and (4.6) into Eqs. (3.13) and (3.14), and equate 
coefficients of powers of a. It follows that the constant \ has the form 


A=Atart+°::,; (4.7) 
where A» , A, , -** are constants. Also, if we recall (2.9) and result (a) in Sec. 2, we see 


that the pressure may be written 


P= : [ Oo dIIg + kXoh + api + -°:, (4.8) 
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* only. In addition, using (4.7) and (4.8), we see that the 
13) yield the equation 


where p, is a function of £, ¢ 
terms independent of a in (3. 


(a + 
24 (Oohys) + == (Oohz)? + Xo = 0, (4.9) 
(0g Of ) 


\ 


whilst the terms independent of @ in (3.14) now cancel. If we consider terms in a in 
(3.13) and (3.14) we obtain 


' | 
hi + 2 fOogee + O:hes + BWkOgh:pre¢0} 


or 
9) 
+2 rer {Oogr + Ohy — WhkOchzeesr} = At p(hypre — Nyogy), (4.10) 
O¢ 
aD, 9d P 
ee + 2 aT; {[hOv(hegre + hogs) + (kOr + Wo)hehge} 
§ $ 
A ; 2 
+9 = [1Ooprere + WOohzegre + (kO, + Wi)hi.} = 0. (4.11) 
df 


We continue discussion of these equations for the special fluid which is defined by 


© = 2p, We=c+II, k= 0, (4.12) 
where yu and c constants, so that 
Oo = 2z, W=ct+/l, k =0, 6, = 0. (4.13) 
Equations (4.9), (4.11) and (4.10) respectively then reduce to 
Suhre + Xo = 0, (4.14) 
Zn + Aipgrrere + 8 = (hzhz-) + 8 - (hyhz.) + 2e{ 2, (hyhye) + = az} = 0, (4.15) 
Sudrre +A, = 4ip(hpyere — hyee;). (4.16) 


If we differentiate Eq. (4.15) with respect to ¢, subtract the result from its complex 
conjugate, and use (4.14), we obtain the following equation for ¢: 


os ; 0° $ 
UUPrrrere + ag? (h-h;«) a ac® (hyshz) = (), 
or 
UWP rrrers + Ghepe(heehey —_ hzhyere) 
oh heehes: = hzhzepere 
4 Shehye(hesheres scone hghgrers) = 0. (4.17) 


5. Flow through elliptical pipe. In this section we consider the solution of Eqs. (4.14), 
(4.16) and (4.17) for the special case of flow through a pipe whose cross-section is the 


ellipse 


~~ 


2 
x 


. } (5.1) 


2 

o Is, 
TS) 
Il 
— 








1956] STEADY FLOW OF NON-NEWTONIAN FLUIDS 307 


In complex coordinates ¢, ¢* this equation may be re-written in the form 
(a? + b*)¢e* — 2a? — B*)(¢? + *’) — 20d’ = O. (5.2) 

At this boundary we assume that the velocity of the fluid is zero so that 
h=0, ¢; = 0, g = 0, (5.3) 


when ¢, ¢* are related by Eq. (5.2). 
Since h vanishes on the boundary (5.2) we try a solution of (4.14) of the form 


h = A{(a’ + b*)¢e* — Ha’? — B*)(¢? + &*’) — 20°b’}. (5.4) 
This satisfies Eq. (4.14) if 
Su(a? + b°)A = —r. (5.5) 
If we substitute h from (5.4) in the right hand side of (4.17) we have 
iu@rrrere — 24A‘a7b7(a* — b*)(¢? — ¢**) = 0. (5.6) 
We require a solution of (5.6) such that ¢, (and therefore g;-) vanishes on the boundary 
(5.2). We therefore assume that 
ing = B{(a? + BW) Ee* — Ha? — B)(g? + f*) — 20°B"}7(¢? — f*’). (5.7) 
This satisfies the boundary conditions and will also satisfy the differential equation 
(5.6) if 





Bi4(a? + 0’)? + (a? — b’)?} — 4A‘07b?(a* — b*) = 0. (5.8) 
Using (5.8), and expressing ¢ in terms of the coordinates xz, y, we have 
644 ‘a7b?(a* — b*)(b?2”? + a’y’? — a’b’)?xy ‘ 
= 2 2\2 2 2\2 . (5.9) 
u{d(a + b) + (a — BY} 


This represents a flow in the z, y plane whose streamlines are the curves ¢ = constant. 
The flow is illustrated in Fig. 1. 











Fic. 1. Secondary flow in a tube with elliptical cross-section (@ > 0). 
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To complete the solution of the problem as far as terms of order a are concerned we 
must now solve Eq. (4.16). In this equation we may put A, = 0 since this only contributes 
terms of the same type as those in the original flow down the tube. The remaining terms, 
which arise solely from the inertia of the material, are found from (4.16), where h, ¢ are 
given by (5.4) and (5.7) respectively. Thus 


8pAB ,, » , . —e 2 
Sugrs = Bp A8 1(@ + OI? — He — CME + 8") 
be 
— 2a*b’}*{2(a? — b*)¢e* 
— (a + b*)(¢? + ¢**)}, (5.10) 
together with the boundary condition 
g=0, (5.11) 


on the boundary (5.2). The solution of this boundary value problem is straightforward, 
but we omit details since the resulting flow will only modify slightly the original flow 
down the tube. 
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—NOTES— 


ON SOME CHARACTERISTICS OF THE FREQUENCY EQUATION OF SMALL 
VIBRATIONS OF HOLONOMIC CONSERVATIVE SYSTEMS WITH STATIC 
COUPLINGS* 

By DANILO P. RASKOVIC (University of Belgrade, Yugoslavia) 

The frequency equation of a holonomic system which performs vibrations about its 

configuration of stable equilibrium is 
A,(A) = |C — AA | = 0, (1) 

where A = (a;,) the square inertia matrix, C = (c,,) the stiffness matrix, \ = w the 
characteristic number. 

When the system contains only static coupling, a,, = 0 for i ¥ k, and in the special 
case when the coefficients of matrices are a;; = 4, ¢;; c, Eq. (1) becomes 


A,(A) = |AI— P| =0, (2) 


where P = A™~'C is the dynamic matrix. 

In the three characteristic cases of the torsional vibrations of light shafts with several 
disks, (Fig. 1), or of the linear oscillations of several coupled masses, (Fig. 2), the matrices 
P have the forms: 











Pp —-? 0 0 
—p 2p 0 0 
. 0 -p 
2p —p 
0 5 =? P) (3) 
2p —p 0 0 | 2p —p 0 0 
—p 2p 0 0 | |—p 2p 
P, = | P=] 
2p —| | 2p —p 
0 0 - —p 2p) | 0 O —p p 
and the polynomials in \ are 
f+) =" + Biv" + +++ + BA + Bo = 0. (4) 


As det P = 0 the matrix P is singular, hence \ = 0 is not a root of the f(A) = 0. 
This condition shows that the generalized momentum of the system is constant; hence 








*Received Jan. 5, 1955; revised manuscript received June 21, 1955. 
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it has only n — 1 degrees of freedom. In the other two cases P is not singular and the | 


systems have n degrees of freedom of vibration. 
The Eq. (2) for these cases can be obtained by recurrence formulas: 


A, = (A — 2p)A,-1 — p’A,-2 = 0 (5) 


and the determinants yield: 


AM = A.) = 0; As? = A + pas) = AS) + pA’ = 0; 6) 
A.” = A,” + pA,” = 0. 

We therefore conclude that clamping of both ends (case b) replaces one disk (mass) 
of the first case; which is the basic one, since Eq. (4) of the other two cases may be 
deduced from the equation of the first case. 

The coefficients B, of (4) are connected by the relations 


B® = BO? — wR” — PB”, + =0,1,--- 


They form a diagonal series of numbers whose differences are A’ = (—1)'(2p)’, A’** = 0, 
where r is the index of the row. In the first case (a) the index r should be replaced by 
r — 1. They also satisfy Newton’s interpolation formula (for h = 1, % = r + 1 or 
2% = r) and may be computed by the formula: 


B.=(1+ A)" = Bo? + ‘ : ")aBs” I scien ‘a _ / + Jara. (7) 


From (7) we obtain the polynomials (4) in the following forms: 


ale 2n — 2 =e 
a) eee | i DN" + ee 


2n — (r ry n—-(r+ n- n> 
a (™ re pe 4 wee + (—1)""'np"™" = 0, 


r 


9 om 
c) ig = . : ‘pa + a 


an — 
+ (" : "px + +++ +(—1)"p" = 0. 
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Putting in the first equation, instead of the index n the index n + 1, one obtains the 
formula for the second case, [1]. 

Contrary to the method of finite-difference equations of second order, [2], the char- 
acteristic numbers in all three cases can be obtained this way. From Eq. (2), by the 
substitution \ — 2p = & = (p*/n) + n, we obtain A, — 7A,_, = (p’/n)"~’ because 
A, = & Ay = 1. Since 7-"A, — 7'~"A,_, = @", @ = p/n, summing from n to 1, while 
A, = 0, one obtains: 

n"A,-1l=qq"-—D/q-1), 1=pexppr/n+t dD), 


and the characteristic numbers are: 


2 pwr : 2 VT : 2 (2v —_ 1)r 
\, = 4p cos on? , = 4p cos an +1)? \, = 4p cos 2(n + 1) (9) 
y=1,2,-:-,n-—1 y=1,2,-°> ,n y=1,2,--: ,n 


Ay > Ag> oo DA; 


Taking into consideration the relations between the roots (9) and the coefficients 
of polynomials (4) we obtain several trigonometric formulas (Table 1). 


TABLE 1 




















Case a b c 
Valeo et — m=n-—!1 ———— * m=n | ( — De —De ;m=n 
x — | fee | 2(2n + 1)’ 
; ’ i ; ~ | 
> CO" gin? « 1 ‘ —(r+ 0) 1 ‘- —(r— ») ae - ") 
— ,silm 2 4 , | 4’ r 4" r 
I] sal 1 ( n ) = - P + ') nti | £ () me 
lll 4-*\n-1/ 4°" | 4&\ n ae 4” \n 4" 
= 1 
» cos 2x 0 0 5 





The symbol =C™ represents the sum of the combinations of the rth class with m 


members, whereas x are the special values of the number 7. 
These relations can be mathematically proved by means of the gamma-functions 


and the geometric series. 
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MOTION OF A VISCOUS FLUID IN A TUBE WHICH IS 
SUBJECTED TO A SERIES OF PULSES* 


By G. N. LANCE (Department of Mathematics, University of Southampton, England) 


Introduction. Consider a circular pipe which contains a viscous fluid. The flow 
due to a constant pressure gradient in a stationary pipe is the well-known Poiseuille 
flow. Now suppose that a pipe, containing a viscous fluid, is subjected to a series of 
pulses in an axial direction which do not appreciably move the pipe. The fluid inside 
will acquire a velocity in the direction of the pulses, relative to the pipe, due to the trans- 
fer of momentum from the boundary by viscous stresses. We determine the flow, relative 
to the pipe, which is obtained when a pressure gradient and a series of pulses act. Since 
the equation of motion is linear the two effects may be considered separately and the 
two results added together. 

In particular we are interested in the case when the pressure gradient and the pulses 
have opposite senses. In this case there arises the possibility that the discharge from 
the pipe could be completely stopped.** In Sec. 1 we consider the simpler case of flow 
in a two dimensional channei. The case of a pipe with circular cross-section is taken up 
in Sec. 3. 

1. Flow in a straight walled two-dimensional channel. Let x be the coordinate 
measured along the channel length, 2a be the distance between the walls and the origin 
of coordinates be midway between the walls. 

For the steady state motion the velocity in the axial direction u(y) satisfies the diff- 


erential equation 


du _— dp 
—= ae —iee (1 
“ dy dz } 
If the pressure gradient dp/dr = —pP (a constant) then the velocity is given by the 
expression 
Pa’ | (4) | 
uly) = a a ) (2) 
Y 2v a 
where the viscous boundary condition of zero slip on the planes y = +a has been used. 


The result may be rewritten 


wn = [1 (Y)] @ 
~ ( 


where u,, = Pa’/3v is the mean value of u across the width of the channel. 


The time dependent equation is 


my —_ 7 pP = 0. (4) 
OY Ot 
*Received March 24, 1955; revised manuscript received August 16, 1955. This paper was written 
when the author was in the Department of Engineering, University of California, Los Angeles, and he 
wishes to acknowledge receipt of a maintenance grant of the F. O. A., Scientific Research Project 
TA01-101-3006 (OEEC 151). 
**This problem was suggested by reports that the engines of certain jet aircraft fail when the guns 


are fired. The failure could be due to fuel starvation. 
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The flow velocity u(y, t) due to the pressure gradient and the action of pulses can be 
considered (since the equation is linear) to be composed of the flow (2) and an additional 
velocity u,(y, t). Introducing non-dimensional variables v, — and 7 defined by 


u(y, t) = uv, T), (5) 
y = ag, (6) 
and 
t=a°T'/v, (7) 
the equation for v(t, 7) is 
2 
ot (8) 


2. Boundary conditions. The boundary and initial conditions which must be 
imposed on v are 


(i) v(t, T) = 0 for T <0 (9) 
N 
(ii) w(+1,7) = >> A,(T —nD) for T>O0. (10) 


In (ii) A, is a non-dimensional quantity giving the strength of the nth pulse and D is 
the time interval between pulses. The differential equation (8) can be solved, subject 
to the conditions (9) and (10), by the Laplace transform technique. If capital letters 
denote the Laplace transform, with respect to 7, of the corresponding small letters 
then we get 





d’v . 
and 
N 
V(+1,p) = >) Ae”. (12) 
Thus 
; cosh tp”? = mile 
V(é, p) = a5 A,e ; (13) 


cosh p _ 


The transform of the quantity of fluid flowing (relative to the channel) per unit 
width of channel, per unit time, is 


a+1 N 
Q(p) = | V(é, p) dt = 2p" tanh p'” >> Aye”. (14) 
J] 0 


The inverse transform of (14) can be found in many tables. (See for example B. van 
der Pol and H. Bremmer [1].) We find 


q(T) = 2 >> A,0.(0, [T — nD]/mH({(T — nD]/n), (15) 


where H(z) is the unit function. 
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The theta function which enters may be expressed in the form 


(0, = 2 > exp | - x(m + ys] 


so that 
q(T) = 4 = A,H((T — nD]/x) > exp } -(m + Ne _ nb) |. (16) 


m=O 


3. Flow in a pipe of circular cross section. We turn now to the case of axially 
symmetric flow in a pipe of circular cross section. The radial coordinate is r and the 
axial direction is x as before. The treatment follows exactly that used in Secs. 1 and 2 
In the present case we require a solution of the differential equation 


dv , lav ow (17) 
a ' £0 aT’ 


where r = aé, ais the radius of the pipe and u,, = Pa*/8v in the present case. The bound- 


ary conditions are 


(i) af, 7) =) fo T <0, (18) 
(ii) v(1,7) = >> A,a(T — nD), (19) 


0 


with the further condition that v is finite on the axis. 
The Laplace transform V(é, p) of v(é, 7’) satisfies 


ev 4 ae pV (20) 


subject to 


V(1,p) = >, Aye”. (21) 


Hence 

, . (ep *) —pnD ‘ 

VE,p) = S45 A,e (22) 
$i 1 “a ip?) >? 
The transform of the quantity of fluid flowing (relative to the channel) per unit time is 
Q(p) = 2m [ eV(c, p) dt = sy (p Ag. (23) 
a ] ee “a n€ a) 
7 J $ o,f $ ip’ (ip 2) 


The inverse transform of Q(p) is 


9 {x —_ ‘ 
(T) = 2x = 45> (3 — 2 > a;” exp [-a(T — nD)\)H(P _ nb)} 
0 m=1 


aT | 2 
or 
N ( ee) 
q(T) = 2x >> A,)2 > exp [—a2(T — nD)]H(T — nD) 
= (24) 


+[3-: 2 > a, exp {—a,(T — nb)} fact - np)}, 


where a,, is the mth root of the Bessel function Jo(z). 
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4. Discussion of the results. In the case of the two dimensional channel the basic 
Poiseuille flow causes a discharge obtained from (3) of amount 2au,,/unit width of channel. 
Thus if we confine ourselves to one pulse the total discharge at time T' is 


ag)? + 4A,H(T/z) >» exp | —(m ot Yr}. 


m=0 
Since the pulses act in a direction opposing the free streaming A, is negative and if T 
is such that 


3 < —A, >> exp | -(m + Vr | (25) 


m=0 
then the flow will be reversed. This inequality can clearly be satisfied if 7’ is small. The 
length of time for which the flow is arrested will depend on the magnitude of the impulse 
coefficient A, . Alternatively, if another impulse follows the first rapidly enough, the 
flow will remain arrested for a further period. 
In the case of the circular pipe the basic flow causes a discharge ra*u,, . So the total 
flux due to a single pulse and the basic flow is 


( be r) 
aval x 4+. 2r Ao)? Z exp (—a,T)H(T) + E —2 b » a, exp (-atr) Jacr}). 
m=1 ad m=1 
It follows, as before that when T is small but not zero the flow will be arrested if 


] — i 
4 < ~Ao >. exp (—a27). (26) 
m=1 
Conclusions. The above analysis shows that it is entirely possible for the flow 
of a viscous fluid, in a pipe or channel, to be arrested when the pipe is subjected to 
one or more pulses of sufficient strength. The time for which stoppage occurs depends, 
essentially, on the magnitude of the pulses, the viscosity of the fluid and the radius of 


the pipe as well as on the frequency of the pulses. 
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AN EXTENDED UPPER BOUND THEOREM ON THE ULTIMATE LOADS 
OF BUCKLED REDUNDANT TRUSSES’ 


By E. F. MASUR (University of Michigan) 


Introduction. This note is intended to supplement an earlier paper® in which the 
collapse load of a buckled redundant rigid-jointed truss was bracketed between a lower 
and an upper bound. As was pointed out at that time, the lower bound theorem permits 
the computation of a sequence of “‘statically admissible’ parameters which approach the 





‘Received August 22, 1955. This paper was written under the sponsorship of the Office of Ordnance 
Research, Contract No. DA-11-022-ORD-1156 while the writer was at Illinois Institute of Technology. 

2E. F. Masur, Lower and upper bounds to the ultimate loads of buckled redundant trusses, Quart. Appl. 
Math. 11,- 385-392 (1954). 
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exact value of the ultimate load within arbitrarily close limits. The upper bound theorem, 
on the other hand, restricts itself to the establishment of a finite (usually small) set of 
load parameters, of which each is known to be too large. The gap between lower and 
upper bound can therefore not be narrowed down arbitrarily; furthermore, the lowest 
of all upper bounds is usually not close enough to the exact collapse load to be of technical 
interest. 

It is easy to see why this should be so. In fact, the proof of the upper bound theorem 
is predicated on the construction of a series of collapse modes which are physically 
unrealistic since they involve the concept of some bars buckling as fixed-end columns, 
while the remaining bars are rigid. This set of collapse modes without joint rotations 
does not include the actual collapse mode; this accounts for the shortcomings of the 
upper bound principle described above. 

An extended upper bound theorem. In what follows, the concept of a collapse mech- 
anism will be extended in such a way as to encompass an infinity of collapse modes. 
With each mode there is associated a “kinematically admissible’ load parameter. 
Moreover, the actual collapse mode belongs to the set of modes to be defined; hence, 
the ultimate load parameter is also kinematically admissible. Finally, it will be the 
object of this note to establish the following. 


THEOREM. The ultimate load parameter is the smallest of all kinematically admissible 
load parameters. 

To prove this principle, let a system of bar forces S, be in equilibrium with the 
external loads identified by the multiplier \4’, and let a set of geometrically consistent 
deflection functions y,(x) exist which satisfy the equations of equilibrium 


(Elyt’)” — Sy’ = 0. (1) 


Let it be assumed further that the change-in-length terms 
aLet 
| s\2 9 
b = 3 (yi) dx (2) 
0 
are derivable from a geometrically consistent set of joint displacements; this is satisfied 
if, and only if, 


> Sh =O & =1,2,---,m). (3) 
; 


Let finally the potential energy associated with the configuration y,(x) be non-positive, 
that is: 


> d a;;0:6; < 0. (4) 
. 3 
When these conditions are complied with, the multiplier 44’ will be designated as 
being “kinematically admissible.” It is readily apparent that the collapse load para- 
meter A> is itself kinematically admissible since, under collapse conditions, Eqs. (1), 
(2), and (3) as well as the equality sign in (4) are satisfied. 
In a manner similar to a previous analogous development’ the satisfaction of the 


3E. F. Masur, A large deflection theory of redundant structures, Proc. Second U. S. Nat. Congr. Appl. 
Mech., 427-436 (1955). 
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first three equations can be shown to be tantamount to the single set of equations 
> Dd a:;.-0,0; = 0 (ry = 1,2, --+ , m). (5) 
7 4 


Hence the load parameter )j’ fulfills the condition of kinematic admissibility if there 
exists a system of internal forces S, (in equilibrium with the external loads associated 
with \4’) and a set of joint rotations 6; such that the relations (4) and (5) are satisfied.‘ 

If now the truss, with its force system so defined, is deformed into a kinematically 
admissible configuration satisfying Eqs. (1), (2), and (3), it follows from the force equa- 
tions of equilibrium at the joints and the compatibility of the bar shortenings 6, that 


w’W = —>> Sid (6) 
k 


in which \4’W represents the work done by the external loads. Similarly, for the same 
configuration, the actual collapse force system S; (corresponding to the collapse para- 
meter A) gives rise to the work relationship 


WwW = —~> Ss, . (7) 
k 


Hence, by combining Eqs. (6) and (7), 


> (Si — S,)d, = (As’ — Aa)W. (8) 


It can readily be demonstrated by means of a few partial integrations and in view 
of the geometric boundary conditions that 


Us + > im S,6, = 4 bB Dd a, 0:8; ’ (9) 
k ‘ ] 


with U, defined as the bending strain energy. On the other hand, it follows from the 
neutrality of the equilibrium associated with S; and Aj that 


Us + >, Sid = 0; (10) 
k 


in (10), the equality sign applies non-trivially only to the actual collapse mode. There- 
fore, upon substitution of (9) and (10) in Eq. (8), 


(As’ — AW = —$ DY De a,;6,0; = 0, (11) 


where the second inequality is the result of (4). 

No generality is sacrificed by postulating a positive parameter \4’; this implies also 
W > 0 by Eqs. (6) and (9), the inequality (4), and the positive definiteness of U, . It 
follows therefore from (11) that 


‘7 
0 


0 2 Xo (12) 


which completes the proof of the theorem. 

In closing it may be pointed out that once an upper bound 4’ to the collapse load 
is obtained, a lower bound may be found, for instance, by reducing all bar forces in 
the same proportion until, corresponding to, say, As , neutral equilibrium is reached. 
It is known then that the actual collapse multiplier \¢ is bounded between Xj and Aj’ . 


‘It is of course to be remembered that the stiffness matrix [a;;] is a function of the force system S, . 
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OVERDETERMINATION OF THE SPEED IN RECTILINEAR 
MOTION OF NON-NEWTONIAN FLUIDS* 


By J. L. ERICKSEN (Naval Research Laboratory) 


1. Introduction. The purpose of this paper is to point out that, in the case of 
rectilinear motion, Rivlin’s theory of non-Newtonian fluids [1] yields two partial 
differential equations for determining a single unknown, the speed. We show that, unless 
the two scalars which occur in Rivlin’s theory satisfy certain definite relations, these 
two equations are independent. When they are independent, it is reasonable to expect 
only a few types of rectilinear motion to be possible. We conjecture that these ex- 
ceptional motions are either rigid or such that the curves of constant speed are circles 
or straight lines. Stone [2] has obtained results in agreement with this for a class of ideal 
materials included in Rivlin’s theory, these materials being very similar to those described 
by the theory of perfectly plastic solids. 

In the case of a fluid flowing through an infinitely long cylindrical tube, to which it 
adheres, it would be natural to suppose that each particle could move with constant 
speed in a straight line parallel to the generators of the cylinder. For this to be possible, 
the cylinder must be a surface of constant speed for some common solution of the two 
equations referred to above. Suppose these are independent. Then, if our conjecture is 
correct, such a flow is impossible except perhaps when the cylinder is made up of portions 
of planes and right circular cylinders. In any event, one may reasonably expect such 
flows not to be possible for cylinders of arbitrary shape. 

Our analysis sheds no light on the question of what type of flow replaces rectilinear 
motion in these tube problems when the latter is mathematically impossible and we 
have found no satisfactory answer to this question in the experimental literature. We 
are of the opinion that it need not be one of the familiar types of flow such as plug flow 
or turbulent flow. 

It should be noted that overdetermination similar to that encountered here arises 
in other nonlinear theories, such as those considered in [3] and [4]. 

2. Equations of rectilinear motion. In Cartesian tensor notation, Rivlin’s equations 
may be written 


tiz.g + pgi = p(dv,/dt + v;,,v,), (1) 
t;; = —pb;; + Fie , Is)di; + F22 , Is)didi; , (2) 
v,,, = 0, (3) 


where #,;; is the stress tensor, g; the body force per unit mass, v; the velocity vector, 
d,; = 1/2(v,;,; + v;,,) the rate of deformation tensor, 27, = — d,;d;;, I; = det d,; , pis 


an arbitrary hydrostatic pressure, and p is the density, assumed constant. The $’s are 
essentially arbitrary functions of the indicated arguments, to be determined by ex- 
periment. 

If each particle moves with constant speed in a straight line parallel to the 2,-axis, 
then v, = v, = 0, dv3/dt + v3.37; = O. If we further require (3) to hold, we obtain 
= 0v;/dt = 0, from which v; = 2f(x, , x2), the factor 2 being introduced to simplify 


V3.3 





*Received Sept. 12, 1955. 
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expressions to follow. Then J, = — f.;f.;, I; = 0, and, from (2), 
Jo 9 Sal fi fafa 0 
twll=—plldell+s%]0 0 S2f+efss. t% 0 f, @ 
Ifa fa 0 0 0 fuhus 


where S; = §,(J, , 0). We assume the body force conservative, g; = ¢,; , and henceforth 
restrict indices to values 1 and 2. Using (1) and (4), we then obtain 


P= (Gof if i): ) P 3 = Gif .3).: ’ (5) 


where P = p — pg is an arbitrary hydrostatic pressure. From (5), P,:3 = P33 = 0, so 
P = ax, + Q(x, , X2), where a is a constant. Then 


Q.; -_ (Gof if i). = Gof i) if. a Gal, ; ’ (6) 


so that the gradient of R = Q + 3/G.d1/, is parallel to that of f. Hence R = R(f) and 


we have 


dR ‘df - (Gof i). = Gale. f 3 + Gof ii ’ (7) 
= (Gif i). _ Gils. if .i 2 Gil ii ’ (8) 


where primes denote differentiation with respect to J, . Equation (7) is equivalent to the 
integrability conditions Q.,. = Q,., and either of these, with (8), gives two differential 
equations for determining f. 


3. Dependency conditions. Suppose that, in some open interval g of values of J, , 
we have 
(GiS2 — GiG2Gi + 21.G1) ¥ 0. (9) 
Suppose further that G, is an analytic function of J, in 9. Under these conditions, we 
will show that, for every a, there exist solutions of (8) such that (7) is not satisfied for 
any choice of R(f). To establish this, we consider a Cauchy problem for (8). Let C be 
an analytic curve, given parametrically in terms of its are length s by 2; = 2,(s), so 
that dx;,/ds dx;/ds = 1. This curve may be chosen arbitrarily, subject to the restriction 
that its curvature not be constant. Let J} < 0 be any value of J, ing and let 
b = (— Ii)’. A Cauchy problem is set by specifying values of f and f.; on C such that 
the compatibility condition df/ds = f.; dx;/ds is satisfied. In particular, we may set 


f =0, f. = bdz,/ds, f= —bdr,/ds (10) 
on C. Differentiating with respect to s, we obtain, on C, 

f 1; dx;/ds = df.,/ds = b d*z,/ds’ = bK dz,/ds, 

f.2; dx;/ds = df ,/ds = —bd’x,/ds’ = bK dx,/ds, 


(11) 


where K is the curvature of C, given by K* = d’x,/ds’ d’x,/ds’. Now (8) and (11) provide 
three linear equations for determining the three unknowns f ,;; as functions of s. Using 
the facts that J,,; f.; = — 2f.:;f.,f.; and that (10) holds, we find that the determinant 
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A of the coefficients in these equations is given by 
A = (dz,/ds)’|S, — 2f7.S1] + (dx./ds)*|G, — 2f7.G1] + 4(dx,/ds)(dx,/ds) ff Si 


a Qf. . Ghana 
= 6 2b wie 


On C, J, = — b’, so (9) implies that A ~ 0. Hence the f,,; are determined uniquely. 
Similarly, all higher derivatives of f can be determined on C, so we can determine a 
power series expansion of f about a point P on C. The fact that, in some circle of non- 
zero radius about P, this power series converges to a solution of (8) can be inferred from 
known results in the theory of partial differential equations. For example, one can refer 
all equations to an orthogonal curvilinear coordinate system in which C is a coordinate 
curve, then make use of the Cauchy-Kowalewsky existence theorem [5, Sec. 94]. It is 
relatively simple to show that the solution thus obtained cannot satisfy (7). On C, 
f = 0 and 7, is constant, so that dR/df, S, , Si , G2 and S: must be constant. If (7) and 
(8) were both satisfied, we could, since $,S: — G2S! # 0, solve (7) and (8) for f.;; and 
I, ,f.; . These quantities would then be constant on C. On the other hand, on C, 


Teh = — 2S iS S..; —2b’*f ,,(dx,/ds)? — 2b’f .(dx,/ds)’ 
+ 4b*f ,. dx,/ds dx,/ds = —2b’f ;; + 2b’f.;; dx;/ds dx,/ds 
= —2b°f ,; + 2b’ df,,/ds dx;/ds = —2b’f,;; + 2b°K, 


where we have used (10), (11) and dzx,/ds dx;/ds = 1. Thus, if /,,;f.; and f,;; were both 
constant on C, we would have b’K constant, which contradicts the assumptions J; < 0 
and K = constant made earlier. Thus f cannot satisfy (7). 

Excluding the pathological case where G, is not analytic in any interval where (9) 
holds, we have left to consider what occurs if the left member of (9) vanishes for all J, . 
We assume S, + 0 which will be the case if, for example, the inequalities discussed in 
[6] are satisfied. If $,$: = S.Gi, we then have SG, = aG, , where a is a constant, or, in 
terms of previous notation, 


$1, ,0) = aF, (J, , 0). (12) 


Whenever (12) holds, every solution of (8) satisfies (7) with dR/df = aa, so f is not 
overdetermined. It then follows that a function f satisfying the basic equations with 
§, = 0 will satisfy them with F, given by (12) for any a. Oldroyd has discussed several 
methods of obtaining solutions for rectilinear motion of materials for which $, = 0 
[7], [8], [9]. It is perhaps worth mentioning that, if (12) holds with a ¥ 0 for all motions 
such that J; = 0, there will exist motions for which the inequalities discussed in [6] are 
not satisfied. If S, + 27,9{ = 0, S, = k(— I.)~'”’, where k is a non-zero constant. In 
this case, the curves f = constant must be straight lines or circles, as has been shown 


by Stone [2]. 
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VECTOR FIELDS ASSOCIATED WITH PLANE PLASTICITY* 
By D. R. BLAND (King’s College, London) 


Synopsis. A symmetric tensor can be associated with a vector field by means of 
Eq. (1). It is found that, for the stress field in plane plasticity, associated vector fields 
can always be found. For certain vector fields, associated stress fields exist. Examples 
are given. 

Transformation between stress and vector. The most general way in two dimen- 
sions that the symmetric stress tensor o,;; can be associated with a vector field v; , when 
no derivatives or integrals of v; are included, is 


o;; = Mé,;; + Now; , (1) 
where M and N are functions of v? . The stress tensor satisfies o;;.,; = 0. Hence, from (1), 
2M'v,v,.; + 2N'v,,. 00; + Nov;,,v; + Now;,; = 9. (2) 
tesolving the vector equation (2) parallel and normal to »; , 
(2M’ + 2N'vi + N)v,vv;.; + Novvw;.; = 0 (3) 
and ** 
€;0,0,(2M'v;,; + Nv;,;) = 0. (4) 
The invariants of the stress tensor can be found in terms of M, N and 2; , 
I, = 0, = 2M-+ Ny (5) 
and 
I, = 0;;0;; = 2M? + 2MNu;, + N*(r})’. (6) 
The stress tensor for a solid in the plastic state satisfies a criterion of the form 
Sh , 12) = 0. (7) 
Substituting from (5) and (6), 
g(M, N,v) = 0, (8) 
" Ressived Sant. 28, 1955; revised manuscript received November 18, 1955. 


**e;; is the tensor «1 = «2 = 0, a2 = 1, en = —1. 
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where g is a known function of M, N and x; . It follows that, with the stress field for a 
solid in the plastic state, an infinity of vector fields satisfying the Eqs. (3) and (4) can 
be associated. M and N are functions of 1; and satisfy Eq. (8); one of M and N can 
always be chosen arbitrarily. 

With a vector field, a stress field can only be associated if the governing equations 
of the vector field can be expressed in the form of Eqs. (3) and (4). M and N are com- 
pletely determined except for an arbitrary multiplicative constant and an arbitrary 
additive constant in M. This determines the associated stress field and yield criterion 
to within two arbitrary constants. 

Two types of vector field are of particular interest, the non-dilatational and the 


irrotational. For the former v;,; = 0; from (3) 
2M’ + 2Nu,+N = 0. (9) 
The other governing equation for v; is, from (4), 
€,,0,0,{Nv;,; — (N + 2N'vp)v;.;} = 0. (10) 
For the latter v;,; = v;,; ; from (4), 
2M’+N=0 (11) 


and the other governing differential equation is, from (3), 


2N'v00;,; + Nv;,.; = 0. (12) 


Examples of associated fields. (i) The irrotational non-dilatational vector field v,,; = 
v;,, and v; ; = 0. By (9) and (11), N = Band M = — 4Buy; + A where A and B are con- 
stants. The associated stress field is o;; = (— }Bu; + A) 6;; + Bo,v; . The invariant 
I, = 2A. This vector field has an associated stress field whose yield criterion is J, = 
constant. 

(ii) Compressible homenergic homentropic fluid flow. The governing differential 
equations are (pv;), ; = 0, p,; + prjv;,; = 0, 4 v5 + S(1/p)dp = hand p = p(p). »; is 
the fluid velocity, p the pressure and p the density. h is constant. From these equations 


v;.; = v;,; and v,v,;v;,; — (dp/dp) v;,; = 0. On comparison with (12), it is seen that the 
associated stress field exists. N(v;) is found by eliminating p and p from dp/dp = — 
N/2N’, vi = 2(h — J(1/p) dp) and p = p(p). The transformation can be expressed in 
the form ¢;; = C(p — po) 6;; + Cpv,v; , where’C and po are constants. This transfor- 


mation has been considered by Inoue [1] and Hill [2]. 
(iii) Shallow water flow at constant head. The governing differential equations [3] are 


vv;,; = — gz.;, (2;),; = Oand g(z. — z) = 4.0; . v, is velocity and z depth. 2, is constant. 
Whence v;,; = 2;,; and v,2,;v;,; — (g2o — 4 %)v;,; = 0. On comparison with (12) it is 
seen that the associated stress field exists. The transformation is o;; = (A + 4 Bg’z’) 
6;; + Bgzv,v; , where A and B are constants. The associated yield criterion is 


(7, - Vz ™ 24) = Boiai( 1, — 3/3* - 2A). 


(iv) Plane strain of a Mises or Tresca solid. The yield criterion is 27, — Ij = K’, 
constant. From (5) and (6), N = k/v;. The associated non-dilatational vector field will 
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be v,.; = 0 and ¢€;,0;0,(v;,; + v;,; = 0 with the transformation 


‘ 


r 


0:;= (z In vi 4 cba + Son, . 
2 Vk 


(v) Plane stress of a Mises solid. The yield criterion is 37, — Ij = 2K’. Substitution 
from (5) and (6) shows that M and N satisfy 


M*? + MNvu; + N°)? = K’. 
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PHYSICAL INTERPRETATION OF PHYSICAL COMPONENTS 
OF STRESS AND STRAIN* 


By DANIEL FREDERICK (Virginia Polytechnic Institute) 


1. Introduction. The purpose of this paper is to give a physical interpretation to the 
several different sets of physical components of stress and strain which have been 
introduced in the literature by different authors and to present and interpret some new 
sets. The question of the physical interpretation of the physical components of stress 
and strain arises when one expresses the equations of mechanics involving these in a 
general coordinate system through the use of tensor calculus. Since the tensor form of 
the equations can be obtained readily, it is of practical interest to have a direct physical 
interpretation of the physical components which are related to the tensor quantities. 

The reader is assumed to be familiar with the elements of tensor calculus as presented 
in the book by Synge and Schild [1] and their basic notation will be used here. All quanti- 
ties are referred to the ordinary physical space which is flat and has a positive definite 
metric form. 

2. The physical components of the stress tensor. These will be introduced through 
the equations of equilibrium using the four sets of components of the stress vectors 
which act on the faces of an infinitesimal curvilinear tetrahedron associated with the 
covariant (contravariant) base vectors g; (g'). When these equilibrium equations are 
compared with corresponding tensor equations, four different sets of physical components 
arise for each triad. 

(a) Physical components of a vector. The four sets of physical components associated 
with a vector Z in general coordinates x’ will be divided into two groups. The first are 
called orthogonal components, since they are obtained by projecting the vector ortho- 
gonally onto the base vectors, and will be denoted by the capital letter O to the left of 
the base letter. Those components along the covariant or contravariant base triads 
which add by the vector or parallelogram law to give the original vector are called 
parallelogram components. They are marked with the capital letter P to the left of the 
base letter. Used as subscripts (superscripts) these letters refer to components along 


*Received July 15, 1955; revised manuscript received October 5, 1955. 








324 NOTES [Vol. XIV, No. 3 


the covariant (contravariant) base triad. The relations between the physical com- 
ponents X,,;, and the tensor components X; , X* are shown in Table I. These were given 
by Ricci and Levi-Civita [2] in 1901. For orthogonal coordinates where the metric 


= 


tensor a;; = 0, for i ~ j, the four sets are identically the same. 


TABLE I* 


Covariant base triad Contravariant base triad 
ahin @ 2 "lan 3" 
Xu = Xa)” Xen = Xe") 


Also, the X,;) components are shown in Figure 1 for plane x’ coordinates. Note: 
Unless otherwise signified the range convention is assumed to operate for all indices. 

(b) Notation for the physical components of stress. Since the stress vectors associated 
with either triad can be resolved into the four components described above, it is possible 


x* curve 





x-curve 


° 
Xw » 


Xo 


Fic. 1. Physical components of a vector in plane 2’ coordinates. 


to introduce eight distinct sets of physical components of stress. The four sets associated 
with the covariant (contravariant) base triad 7',(7”) will be denoted by placing (7j) as 
subscripts (superscripts); while the letters O and P to the left of the base letter T have 
the meanings described above. 

(c) The pT.;, components. This set of physical components arises upon resolving the 
stress vectors associated with the triad 7’, into parallelogram components along the 
covariant base directions. In order to relate these to a stress tensor, the equation of 
force equilibrium along the z* curve will be written and compared with a corresponding 
tensor equation. 

To carry this out, consider the infinitesimal tetrahedron associated with the triad 
T; shown in Figure 2, where the face perpendicular to the unit normal vector 7 will 





*The summation convention is not used here. 
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be dS and the face normal to g‘ will be dS,,;,. The stress vector on the face dS[dS,,,] is 
i{Z,,,]. These elemental areas are related by the formula 


dS,;, = n{a‘’)'dS. (1) 


Also, consider the components of force p7';;;, dS,;, on the face dS;;). These three 
forces, for 7 = 1, 2, 3, are projected onto the 2* curve by multiplying by the cosine of 


x curve 


ta 








t 


a 







~ x curve 


x'- curve 


tay 


Fic. 2. The curvilinear tetrahedron with faces normal to the contravariant base vectors. 


the angle between the x’ and x‘ directions which is a,; (a,,a;;)~'’* and summing over j. 
The result is 


3 
yo Puttin “ (2) 


7=1 
Summing over 7 in Eq. (2) gives the sum of the force components on all three faces 
dS,,, in the direction of the x* curve and this must be equal to the component of force 


on the face dS in the 2“ direction, which is of,,, dS. Hence 


3 
pT (5j)08 50x (Qua;;)) = otaydS. (3) 
t,7=1 
Using Table I and Eq. (1), this becomes 


3 


t, = > FP td hi. (4) 


i,j=1 
Comparing this with the tensor equation 
t, = Tn; (5) 


yields a relation between the mixed tensor components and the p7';;, physical com- 


ponents 
3 
T'. oe + # pT ¢:;)Qx;(a"*/a;;)' a (6) 
i=l 


By raising the index k, there results 


T= pT i;(a"'/a;,)"”. (7) 
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This is tabulated in Table II with the seven other relations which can be found by 
a similar procedure. The numbers refer to the reference where the relation was given. 


TABLE II 


Covariant base triad Contravariant base triad 
T = pT ,;(a"'*/a;;)'  [3][4] TT’, = pT"(G,;/a;)"" [5] 
YY «= Faia P , ea 
T= "Fade" Tt. = Te fa 
Mae Teo .6|CUG = Te Ga” 


For the special case of orthogonal coordinates these eight sets of physical components 
become identical. 

3. Physical components of infinitesimal strain. The material of this paragraph is 
presented in a form which presumes that the reader is familiar with the notation and 
presentation of strain as given by Green and Zerna [4]. Assuming infinitesimal strains, 
these authors present the physical components 


Vii , 
4%) — (8) 
YX Qs 
and, for i ¥ j, 
TH = [27:5 —" 0: (¥ ci) + Voi) |[as.a;; _ (a;;)"] ee (9) 


In these equations y;; is the covariant strain tensor and a,; is the metric for the unde- 
formed state. The physical components y,;;, are changes of length per unit of length for 
line elements along the coordinate curves, while y,;;, is the change of angle between the 
covariant base vectors g; and @; . 

Another useful set of physical components of infinitesimal strain can be obtained by 
considering the line elements along and the angles between the contravariant base 
vectors. Following the procedure used by Green and Zerna for y;;;) and ¥;;;) , the final 
expressions for the corresponding strains y‘''’ and y"*”’ are the following: 


(ia) y"" \ 
= 1%, 10 
Y a (10) 
7 -_ [2y" a ie. ai + y")][a‘*a’? a. yr”. (11) 
For orthogonal coordinates y"""’ = ys) and y°" = Yi 


4. Conclusions. For stress (strain), there are eight (two) convenient sets of physical 
components which might be used for solving boundary value problems in elasticity 
involving general coordinates. 

Although the results are presented for a triad tied to the x” coordinates, the same 
ideas can be extended to a general point using the theory of non-holonomic coordinates 
as explained by Truesdell [6]. 
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A NOTE ON THE DECOMPOSITION OF STRESS AND STRAIN TENSORS* 
By BERNARD BUDIANSKY anp CARL E. PEARSON (Harvard University) 
A well-known theorem of vector analysis states that an arbitrary vector field can be 


decomposed into the sum of two fields, one of which is irrotational and the other solenoidal. 
Analogous theorems for stress and strain tensor fields are presented in this note. 
Theorem 1. A symmetrical (stress) tensor o,;; defined in the region V bounded by 
the surface S may be written 
a) ies oi; + ot; ’ 
where o/,; and of; have the following properties: 
(a) Otis ao ai fi in V, 
ain; = T; on a portion Sz of S, 
where f; and 7’ are prescribed’, and n, is the unit normal to S. 


(b) the (strain) tensor e{/ derived from o{} by the Hookean relation 


di = ll) = £0 +e — wl 


es 


is related to some (displacement) vector field u{’ by 
Yt - = , ” 
“ * (us! ; + us’), 
° 0 Y ry 
where u;{’ takes on the prescribed value u; on S, = S — Sz. 


In other words, the theorem states that any stress field can be decomposed into the 
sum of two fields, one of which obeys prescribed conditions on internal and surface equilib- 
rium and the other of which provides (through Hooke’s law) strains that satisfy internal 
and surface compatibility requirements. The analogy with the corresponding result 
for vector fields is evident: the equilibrium condition involves the divergence of the 
stress tensor (a solenoidal vector is divergence-free) and compatibility requires that 





*Received December 20, 1955. 
‘Here f? and 7°f have the character of body force and surface traction, so that if Sg is the entire sur- 


face, they must be so chosen as to satisfy overall static equilibrium. 
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strain be derivable from displacement (an irrotational vector is derivable from a 


potential).* 
The proof of the theorem rests on the fact that o// may be chosen as the solution to 


the following mixed boundary value problem for o/} , ui’ , and €(; : 


oe =cut ff 86 in, 


ff = Lot) = Ful’; + uj’) in V, 
0 nN. O;n; =— T; on Sp 9 
Tad =U; on S4 
Then the choice ¢/; = o;; — o// satisfies 
gh = —f; in V, 
aim, = T; on Sz, 


and so o/; and 4} constitute the elements of the required decomposition. 
Theorem 2. A symmetrical (strain) tensor e;; defined in V may be written 


; Lad 
€5; ed + €5; ’ 


where ¢/; and ¢;; have the properties 


(c) &f = 3(ul’, + ut’;), where u/’ is a (displacement) vector that takes on the values 
u; on S,, 
(d) the (stress) tensor o/; = H(e/;), where H is the inverse of ZL, satisfies 
oi34 = —f; in V, 
oim, = T; on Sz, . 
This theorem follows immediately from the application of Theorem 1 to the stress 
tensor c,; = H(e;;) with the result that the decomposition ¢;; = o/; + o/; obeys con- 
ditions (a) and (b). Then the strain tensors ¢/; = L(o/;) and ¢/} = L(o‘t}) provide the 


elements of the required decomposition of ¢;; . 








The reader is cautioned not to infer that o// satisfies the conventional Beltrami-Michell compatibility 


equations, since these involve equilibrium conditions that o// need not satisfy. 


ON VARIATIONAL PRINCIPLES AND GALERKIN’S PROCEDURE FOR 
NON-LINEAR ELASTICITY* 


By BERNARD BUDIANSKY anp CARL E. PEARSON (Harvard University) 


Variational principles are discussed by Greenberg [1] for the following problem of 
non-linear elasticity** in the region V bounded by the surface S = S,4 + Sz: 


6:3. =0 in V, 

€ 5; = (u, ; + u;.;) in V 
U; =, on S,z , 

o,n, = T; on Sz, 


*Received December 20, 1955. 
**The non-linearity is incorporated only in the stress-strain relations; displacements are still as- 


sumed small. 
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where 
“= €:;(011,012, see), 
> 4; (€11,€12, s+), 
The unambiguous formulation of the variational principles requires that the quantities 
o,,de;; and ¢,,do,;; be exact differentials; it is accordingly assumed that d0,;/d¢,, = 
dc,,/¢€;; and that de;;/dc,, = d¢,,/d0;; . Then it is shown in (1) that 
T. If u; is a solution, with e;; = 3(u;,; + u;,,), then 


| | [ tibinta Vial ~ i te as | =® (1) 
JV JO Sp 


for all du; satisfying du; = 0 on S, . Here the upper limit in the inner integral of the 
first term is the ¢;; state corresponding to the u; solution. 
II. If o;; is a solution, then 


| | [ €:j(011,012, +**) do;; dV -f Uso «iN; as, | =0 (2) 
Jv Jo Sa 


for all 5c,;,; satisfying é0;;,; = Oin V and 6e,;;n; = 0 on Sz . The upper limit of the inner 
integral in the first term is the o,; state of the solution. 

In the practical application of these principles the procedure is to enforce (at least 
approximately) the satisfaction of Eqs. (1) and (2). The question therefore arises whether 


the following inverse theorems are valid: 

III. If (1) holds for a particular u; satisfying u,; = ui on S, , and for all du; satisfying 
du; = 0 on S, , then this u; is a solution. 

IV. If (2) holds for a particular o,; satisfying o;;,; = Oin V, o;;n; = Tj on Sz , and 
for all d0;; satisfying 60,;,, = Oin V, 60;;n; = 0 on Sz, then this o;; is a solution. 

As will be shown, the validity of III can be readily established by elementary means, 
but the proof of IV is less obvious. To prove III, we note that the volume integral in 
(1) may be transformed as follows: 


5 [ — or | eutin 4 « [ o,;8u,,, dV 
Jy Jo v v 
= [ (o.;5u,),,dV — [ o:;,<du; aV. 
Then use of Green’s theorem and substitution back into (1) gives 
-[ ¢<;,48, AV + [ (o,, — T*) bu, dS, = 0. 


It then follows from the freedom of choice available for du; that o,;,, = 0 in V and 
o,;n, = T°? on Sz , and so, since u; = u} on S, , it must be true that u,; is indeed a 
solution. 

The verification of IV does not seem to be possible in so direct a fashion and so 
recourse will be made to the introduction of the following auxiliary theorem: 

V. If ¢;; is a (strain) field in V that satisfies the requirement 


| €;0;; dV — Uso 4M, dS, = 0 (3) 
Vv Sa 
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for all o;; that obey the admissibility conditions o;;,; = 0 in V, o,,n; = 0 on Sz , then 
¢;; is derivable from a displacement field u; that equals u on S, . 

The proof of V depends on the proposition proved in [2] that e,; can be decomposed 
into €/; + ¢/} , where ¢/} is derivable from a u{’ that equals u; on S, , and ¢/; provides, 


through Hooke’s law, a stress o{; = H(e!,;) that satisfies o/;,,; = 0 in V and a/;n; = 0 


on S; . Then (3) may be written as 


| éf,0;;€aV + [ éf'a;; dV — [ uyo;n; dS, = 0. 
JV S 


JV SA 
By the principle of virtual work the last two integrals cancel for all admissible o;; , and 
the admissible choice of o;; = H(e{;) then gives 


[ este.) av = 0. 
- 
But since the integrand is necessarily positive for e{; ~ 0 (since it is essentially the 
strain-energy-density of linear elasticity) it follows that e{; = 0. Hence ¢e;; = €// , and 
so has the properties stated in V. 

The proof of IV is now immediate; from (2) 


| €;;60;; dV — | u;5o;n; dS, = 0 


Li JSA 


for all admissible 5c,;; . Hence, by V, ¢;; is derivable from a u, that equals u; on S, , and 


therefore o;; is a solution. 
The auxiliary theorem V also provides a theoretical basis for the use of a Galerkin- 


type procedure in the approximate solution of non-linear elasticity problems. Such a 
‘procedure could consist of the following steps: 
(a) Assume an approximate solution of the form 
N 
of = 64; + > —_ 


where 


)) ry0 y 
o3,n, = T;, on Sp 

k) . . , P * 
oi;n; = 0, on Sz fork = 1,2,3,---,N 


and where the constants a, are to be chosen so that 


(b) ff eslohoh, «MP av — [ rn, dS,=0, k=1,2,3,---,N 


SA 


k : . . . (k ° y 
where the \;;' are arbitrary symmetrical tensor fields that satisfy \;;\; = 0 in V and 


Nn, = @QonS;. 

The d; need not coincide with the o;;’; if however they do coincide, this Galerkin 
procedure becomes equivalent to a Rayleigh-Ritz procedure based on (2) and the same 
expansion o# . 

All of the above results may be considered applicable to theories of plasticity of 
the deformation type if no “unloading” takes place (see [1]); in addition the various 
theorems can easily be generalized to the case of non-zero body force. 
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NOTE ON THE EQUATIONS OF SHALLOW ELASTIC SHELLS* 
By P. M. NAGHDI (University of Michigan) 


In the following, a system of differential equations is deduced for thin shallow 
elastic shells (of uniform thickness) with small displacements, which include the effect 
of transverse shear deformation. The corresponding equations of the classical theory, 
where the effect of transverse shear deformation is neglected, are contained in the works 
of Marguerre [1] and Green and Zerna [2], and have also been employed recently by E. 

teissner [3, 4]. Although the equations sought may be obtained (through appropriate 
approximations) from the results given in [5], considerable space will be conserved if 
the basic equations of the theory are referred to Cartesian coordinates. 

Let the equation of the middle surface of the shell be written in the form z = z(z, , 22) 
and let:** w denote the displacement in the z-direction and u,; and 8; be respectively the 
displacements of the middle surface and the change of the slope of the normal to the 
middle surface, along the Cartesian coordinate axes x, and x, . With the notation N,; , 
M,, , and V, , for the stress resultants, stress couples, and transverse stress resultants, 
respectively, the stress differential equations of equilibrium are 


Nis tpi = 9, (la) - 
Vit 2iNaul. +a = 0, (1b) 
Vi; = Mi;.;, (1¢) 


where comma denotes partial differentiation and p, , p. , and q are the components of 
the load intensity in x, , 2, , and z directions, respectively. 

The stress strain relations, which include the effect of transverse shear deformation, 
may be written as 


65 = 3[ui.;5 + uj.) + (z,;w,; + z,,w,;)] 


' (2a) 
= C —vNin bj; oa (1 ao v)N;,;], 
M;; -_ $D[2vB x45: + (1 = v)(B;,; + B;.)], (2b) 
ae ’ _ | ws 
B; 8 5Gh i a) (2c) 


*Received December 5, 1955. The results presented herein were obtained in the course of research 
sponsored by the Office of Naval Research under Contract Nonr-1224 (01), Project NR-064-408, with 
the University of Michigan. 

**Throughout this note, the Latin indices 7, 7, and k have the range 1 and 2 only. Repeated indices 
imply summation over all the values the indices may take. Whenever one of the repeated indices is 
placed in parentheses [as in Eq. (3)], the summation convention is suspended for that index. 
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where ¢;; are the components of strain at the middle surface; 5,;; is the Kronecker delta; 
h denotes the shell thickness and is assumed to be uniform; FZ, G, and v are Young’s 
modulus, shear modulus, and Poisson’s ratio, respectively; C = Eh and D = Eh*/ 
(12(1 — v”)]. 

Introduction of Airy’s stress function F into Eqs. (la) by the relations 


Na = (vr -- | Di) dein )b. —F ;;, (3) 


where V’* = 0°/dx,0xz, , and their substitution into the relevant compatibility equation 
yields 


V'V'F = C422, 1212 Z.11W,22 — 2,228 11} 


bs {L/ Ps a, | sa If P2 az, |. - wish (4) 


which in various forms may be found in [1], [2], and [3]. 
From (2b, c) and (1c), there follows 


1D{(1 — vy) VB; + (1 + »)B;.:;} — 8Gh(B; + w,,) = 0, (5) 


which can also be expressed in terms of V; and w. Next, substitution of (2b) into (1b, ce) 
with the aid of (2c) and (3) results in 


DV‘w — (1 —AVI(V OV’ — 2.4;5F 43] = (1 — AVG 


7 [= [ > ax, | + le |p, az, |} (6) 


where \ = A?/5(1 — v) and V* = YV’V’. It may be noted that elimination of the function 
F between (4) and (6) yields an eight-order partial differential equation in w alone, 
namely 

= eo 3” 
DV‘*7‘*w — (1 - Av)| (W'NV — 25357 


Ox; OX; 


|prien side of Eq. (4)] 
(7) 
= VY*[Right side of Eq. (6)], 


which for cylindrical shells, upon the neglect of transverse shear deformation (A = 0), 
reduces to the corresponding equation given by Donnell [6]. 

Equations (4), (5), and (6) constitute the required differential equations for shallow 
shells. An alternative formulation in terms of displacements may sometimes be pre- 
ferable; this is achieved by reducing the system of Eqs. (1) and (2) to five partial differ- 
ential equations in u, , 8; , and w. Thus, in addition to (5) we have 


ee {(v%« + 7 <a seal + | (07% + 2,,77w) 
a / 


2(1 + »v) 
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and 


Y 


oa + i=» [z, (Vu; 


l-v 2 


$+ 2,,77w) + 2,c5(tis + uy.) +P 270s + i (9) 


DV‘w -— (1 - rv7»/ 





+2,,w,;,+2,.0,,,)+0 —- 2.2.5.5) = [Right side of Eq. (6)]. 


In conclusion, it may be mentioned that as in [5] for vibration problems of shallow 
shells, the effect of rotatory inertia can be easily added to the above equations. 
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ON THE ERRORS IN ANALOGUE SOLUTIONS OF HEAT CONDUCTION 
PROBLEMS* 


By JIM DOUGLAS, JR. (Humble Oil & Refining Co., Houston, Texas) 


1. Introduction. Many practical heat conduction questions lead to problems not 
amendable to the methods of classical mathematical physics. Consequently, numerous 
methods for approximating the solution have been developed; many of these procedures 
result from the replacement of the differential equation by either finite difference equa- 
tions, which are usually solved by means of a digital computer, or by systems of ordinary 
differential equations, which are often treated using an analogue device of some kind. 
Our interest will be confined to the latter of these types. 

Useful analogue machines for the solution of heat conduction problems include 
differential analyzers [3, 4] and resistance-capacitance circuits [1, 6]. In either case, the 
derivatives with respect to one variable are retained and all others are replaced by finite 
differences, and the resulting system of ordinary differential equations is solved. For the 
heat flow equation in one space variable, it is possible to retain either the space derivative 
or the time derivative; however, keeping the time derivative is much more desirable 
from a least work point of view [4]. With more than one space variable, this is the only 
approach that can treat the space variables in a symmetric manner. 

The replacement of the space derivatives by finite differences results in an error in 
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the solution of the original problem. Hartree [4] has conjectured that the magnitude of 
this error is O(Az’), where Az is the greatest space increment. The object here is to give 
a simple proof of this statement. The argument will be an application of a previous result 
on the backwards difference equation [2]. An interesting point is that the difference 
equation result was obtained by adapting a convergence proof given by Rothe [7] for 
the difference-differential equation obtained by replacing the time derivative by a finite 
difference. 
2. Analysis. Consider the boundary value problem 


ou Ou Ou 
—-; — = —— (x, > < mm 
da * ay? at’ z,yeR,O<t<T, 


u(x, y, 0) = f(z, y), 





u(x, y, t) = g(x, y, b), (x,y2eS,0<t<T, 


where R is a bounded, connected set in the plane and S is its boundary. S will be assumed 
to be sufficiently regular in shape that it can be given by a finite number of polygons 
through grid points of a lattice of points of the form (tAz, jAy) for some Az and Ay. 
Let z,;;(#) represent the value of a function of ¢ attached to the point (Az, jAy). Then, the 
analogue solution w;;(¢) satisfies the following system of ordinary differential equations: 


wos as . lw; ;(t) ; ; * 
Arw;,(t) + Ajw,,() = uth" ; (t(Az, jAy)eR,O<t<T, 


w;(t) = g(tAz, jAy, b), (iAr, jAy)e S,O< t <T, 
w;;(0) = f(zAz, jAy), 


where A?w,,(t) is the divided second difference (w,+,,; — 2w,; + w;-1,;)/(Axr)’. Provided 
g(x, y, t) is a continuous function of ¢ for each (x, y) e S, it is well known [5, p. 71] that 
2.2) possesses a unique solution. Moreover, if g, exists and is bounded uniformly on 
(x, y) e Sfor0 <t < T, then w,;(t) possesses a bounded second derivative for0 <i < T. 


In this case, 


i ‘ w,,(t) — w,,(t — Abt) 1 d’w;; 
hie, Ad + Aba, AQ aw SA Ea 5. = Se ae (2.3 
iD HF AD At 2 dt ) 
with d’w,,;/dt? being evaluated at some point in (¢ — At, t). Let (2.1) have a solution 
u such that u,, , Uzz2: , and u,,,, exist and are bounded for (x, y) e RandO < ¢ < T; then 
2 2 u;;(t) — u;,(t — Ab) 1 du 
Avu;;(t) + Avu,;(b) —_—— ——< Ai 
2; (0) yi; At v 95 OF 
(2.4) 
1 o*u . 1 fu ; 
oS tas! as oon \4 
F 79 azt 49 + 79 ay (Ay)’. 
Let 
(2.5) 


V() = u.,(0 — w,,(0 
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denote the truncation error incurred by the use of (2.2). Then, assuming Ay = 0(Az), 


v;;(t) _ v(t — At) 
At 





Aiv; ;(t) + Ai; ;(t) —_ + A,,()At 


+ B,,(t)(Aa)’, (iAz, jAy) eR, 





r (2.6) 
v(t) = 0, (tAz, jAy) e S, 
v;,(0) = 0, (tAx, jAy) e R, J 
where A,;(t) and B,;(t) are bounded. If A = maz | A,;(t) | and B = maz | B;;(é) |, then 
| v(t) | < ATAt + BT(Az)’, (2.7) 
as 
max | v,;(t) | < max | v,,(¢ — Ad | + A(At)? + B(Az)’At (2.8) 


by the growth lemma of [2]. As v,;(¢) is independent of At, let At diminish to zero, holding 
Az and Ay constant. Thus, 


| v(t) | < BT(Az)’, (iAr, jAy)eR,O<t<T. (2.9) 


This completes the proof of Hartree’s conjecture for the cases satisfying the hypotheses 
made above. While for simplicity the proof has been given for the linear equation (2.1), 
it is apparent from [2] that the same argument can be applied to the quasi-linear equation 


Ou au Ou 
ax? sith T 332 _ F(z, , ts ae * u) = + Ga never Me, 48 (2.10) 
if it can be shown that d’w,,....,;,,/dt’ is bounded as needed above. 
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BOOK REVIEWS 


(Continued from p. 276) 


50-100 Binomial tables. By Harry G. Romig. John Wiley & Sons, Inc., New York, and 
Chapman & Hall, Ltd., London, 1953. xxvii + 172 pp. $4.00. 


The value of the individual terms and the sum of the first xz terms of “[P + (1 — P)]” are given 
for n in the range 50 to 100 in steps of 5 and p in the range 0.01 to 0.50 in steps of 0.01 to six decimal 
places although the last place is doubtful. The introduction, written in a clear and concise manner, 
describes the use of the tables, its application to quality control and its relation to the ratio of the in- 
complete 8-function to the complete 6-function. The tabulated data are easy to read but the nature 


of the entries prevents a more uniform tabulation. 
S. L. Levy 


Stress waves in solids. By H. Kolsky. At the Clarendon Press, Oxford, 1953. x + 211 
pp. $5.00. 


This book, the thirteenth in the Monographs on the Physics and Chemistry of Materials series, is 
a welcome addition for those interested in the propagation of stress waves in solids. The first part of 
the book treats the classical theory of stress waves in elastic materials and unites many topics which 
have only appeared in the literature. It is unfortunate that the author did not insert the the results 
for the reflection and refraction of a stress wave at a plane boundary between two anisotropic materials. 

The second portion of the book treats “imperfectly elastic’’ materials and contains a well written 
account of the present state of affairs of internal friction in solids. Chapter VI, Experimental Investiga- 
tion of Dynamic Elastic Properties, presents the various techniques used in the measurement of stress 
waves in solids and discusses the results of the several types of measurements. The last two chapters 
of the book treat plastic and shock waves in materials and the experimental results in this field—a 
large portion contributed by the author. 

It is the reviewer's belief that the author has accomplished his aim stated in the introduction that 
the book should be easily read by anyone trained in mathematical physics. 


S. L. Levy 


French bibliographical digest. Science: Mathematics (Part I: Pure Mathematics). Cul- 
tural Division of the French Embassy, 972 Fifth Avenue, New York 21, N. Y. 
iv + 128 pp. 


This publication contains abstracts or reviews of the most significant contributions made by French 
mathematicians between 1949 and 1954, and aims at making French mathematical work better known in 
the United States. 

The present volume is devoted to pure mathematics; a separate volume on applied mathematics is 
to follow. 

Copies of these publications will be available, on request, without charge to libraries, university 


departments, and individual scientists. 
W. PRAGER 























